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Abstract 

Suppose  we  observe  a  realization  of  size  n  of  a  Gaussian  stationary  sequence  and  we  estimate 
dj,  a  functional  of  the  spectral  density,  by  9X.  For  example,  9j  and  9X  might  be  the  true  and 
the  empirical  rate  at  which  the  process  crosses  a  certain  level.  It  is  natural  to  ask:  (i)  can  the 
sampling  distribution  of  9X  be  consistently  estimated?  and  (ii),  can  a  better  estimator  of  9/ 
than  9X  be  constructed  on  the  basis  of  9X"! 

In  this  paper  we  describe  approaches  which  allow  answering  both  questions  affirmatively  in 
certain  situations.  The  approach  is  based  on  resampling  the  data,  i.e.,  using  the  original  data  to 
specify  a  distribution  from  which  new  samples  are  drawn.  In  the  first  problem  the  data  induce 
an  estimate  of  the  spectral  density  /„  and  the  sampling  distribution  of  9X  under  /„  is  used  to 
estimate  that  under  /.  We  establish  results  for  estimators  of  well-behaved  functions  of  linear 
functionals  of  the  spectral  density.  In  the  second  problem  we  pretend  the  data  are  circular  so 
that  the  periodogram  of  the  data  is  a  sufficient  statistic  and  thus  the  conditional  expectation  of 
9X  is  a  better  estimator.  If  the  data  are  not  really  circular,  this  mechanism  introduces  bias  but 
the  reduction  in  variance  may  be  substantial  enough  to  reduce  the  mean  square  error. 
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Chapter  1 

Introduction 


1.1  Setup  and  problems 

Suppose  we  observe  Xg  ~  N( 0,  £n),  *  realisation  of  sise  n  of  a  Gaussian  stationary  process 
{X,}‘f00.  Since  the  process  is  stationary  £„  is  a  Toepliti  matrix  and  completely  specified  for 
all  n  by  a  spectrum  /,  and  we  may  write  Xg'1  ~  N(0 , /).  Furthermore,  suppose  we  estimate 
dj  —  6(f),  some  functional  of  /,  by  procedure  9Z. 

In  this  paper  we  explore  the  use  of  resampling  methods  in  connection  with  two  problems: 

1.  consistent  estimation  of  the  sampling  distribution  of  9Z  and 

2.  using  9Z  to  find  a  better  estimator  of  9f. 


1.2  General  strategy 

The  basic  strategy  behind  our  approach  to  each  of  these  problems  is  based  on  resampling  the 
observed  data,  i.e.,  using  the  observed  data  Xg~l  to  specify  a  distribution  from  which  we  then 
draw  samples  TJ*-1  which  axe  sometimes  called  bootstrap  samples. 

There  axe  two  steps  implicit  in  this  approach:  the  first  defines  the  procedure  through  which 
the  original  data  specify  the  distribution  from  which  new  samples  are  to  be  drawn;  the  second 
indicates  how  these  bootstrap  samples  are  to  be  used  to  solve  each  of  the  two  problems  at  hand, 
e.g.,  using  Monte-Carlo  to  approximate  distributions  and  their  moments. 

1.2.1  Strategy  for  problem  1 

For  the  first  problem  we  proceed  as  follows: 

1.  Use  Xg~l  to  estimate  the  spectrum  /  by  /„  via  some  procedure 

2.  Draw  Yg~l  ~  N(0,fn)  and  calculate  0,  for  each  sample 

3.  Take  the  Monte-Carlo  approximation  of  the  (conditional  on  X” ~')  sampling  distribution 
of  9y  as  the  estimate  of  the  sampling  distribution  of  0«. 

1.2.2  Strategy  for  problem  2 

For  the  second  problem  we  proceed  as  follows: 

I.  Use  X”~l  to  construct  a  sufficient  statistic  Sn  -  S(Xg~l). 
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2.  Draw  samples  K0n~'"  ~  C(X£  l|S„)  and  calculate  9y  for  each 

3.  Approximate  0*  =  £(0x|SB)  by  the  average  of  the  0y’s. 


1.3  Results  for  problem  1 

In  order  to  be  able  to  give  fairly  exhaustive  coverage  to  problem  1  we  will  restrict  attention  to 
estimates  of  linear  functionals  of  the  spectral  density  and  well-behaved  functions  of  these. 


1.3.1  Linear  functionals 

Linear  functionals  of  a  spectrum  /  are  of  the  form 


9,  =  J  6(*)f(w)dw 


and  they  can  be  efficiently  estimated  (cf.,  e.g.,  Hasminskii  and  Ibragimov,  1986)  by  the  corre¬ 
sponding  linear  functional  of  the  periodogram  J*  of  the  data 


9Z 


/ 


9(w)Il{w)du>. 


For  example,  covariances  and  their  MLE’s  (under  the  Toeplits  model)  are  of  this  form. 


1.3.2  Results  for  linear  functionals 

It  simplifies  matters  to  estimate  the  sampling  distribution  of  Un  =  y/n(9x  -9f)  rather  than  that 
of  0*.  As  mentioned  above,  we  first  use  Xq~'  to  estimate  the  spectrum  /  by  /„  and  then  draw 
samples  Y0n_l  from  £(Y0n_l  |A£~l)  =  JV(0, /„).  If  we  let  Vn  =  s/n(9y  -  9f-J,  we  can  then  use 
C(Vn\X£ -l)  to  estimate  C(Un). 

The  device  used  in  deriving  our  results  is  to  construct  Un  and  Vn  on  the  same  probability 
space  with  the  property  that  C(Un)  =  C(Un)  and  -1)  =  £(V„|Xo~l). 

Our  main  result,  Corollary  3.5,  establishes  that,  under  regularity  conditions  on  9,  /,  and  0n, 

dp(C(Un),C(Vn \xrl))  =  oP(l), 

i.e.,  the  Prohorov  distance  between  the  distribution  of  interest  and  its  estimate  is  small  with 
high  probability. 

The  proof  of  this  result  is  an  easy  consequence  of  Theorem  3.4  which  establishes  Un  -  Vn  - 
o,(l)  under  the  following  regularity  conditions: 

1.  /,  9  satisfy  a  Lipschits  condition. 

2.  11/  -  Ml  =  fo  I/M  -  /nHI’dw  =  0,(1) 

3.  P(infu  /n(u»)  >  0)  -  1 

4-  *«Pw  AM  =  Or(  1) 

5-  supw/-l(<*>)  =  Op(l) 

6-  *upw  |/4(w)|3  =op(n/ log5  n). 
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We  then  extend  these  results  to  differentiable  functions  h  of  linear  functionals  in  Theorem 
3.8  and  Corollary  3.9,  in  which  case  the  variables  of  interest  are 

u:  =  sa  (m«.i  .  ■  •  • .  «.P)  -  m*/i  . 

and 

^  (*(«,. . «W)  ~  *(«/.!.-  •  '  .*/„,))  • 

1.3.3  Spectrum  estimation 

The  motivating  idea  for  assessing  sampling  distributions  is  a  practical  one:  we  would  like  to  use 
this  method  in  real-life  applications,  for  example,  to  estimate  the  sampling  bias  and  variability 
of  statistics  of  interest.  We  mentioned  in  the  previous  section  that  this  resampling  strategy  will 
only  work  if  the  estimate  of  the  spectral  density  used  to  generate  the  bootstrap  samples  satisfies 
certain  properties. 

In  our  simulations  we  implemented  the  spectrum  estimation  technique  derived  in  Wahba 
(1980)  which  estimates  the  logarithm  of  the  spectral  density  by  smoothing  the  log  of  the  peri- 
odogram. 

We  chose  Wahba’s  procedure  for  two  reasons:  (i)  there  is  an  easily  programmable  objective 
criterion  for  choosing  the  amount  of  smoothing  which  makes  the  procedure  automatic  and  thus 
suitable  for  large  simulation  studies;  (ii)  the  degree  of  differentiability  of  the  estimate  can  be 
controlled. 

The  amount  of  smoothing  and  the  degree  of  differentiability  of  Wahba’s  estimate  are  con¬ 
trolled  by  two  tuning  parameters.  Choosing  the  value  of  the  smoothing  parameters  is  essentially 
a  choice  between  controlling  bias  or  controlling  variability  (and  hence  regularity  properties)  of 
the  resulting  estimate  of  the  spectral  density. 

Specifically,  let  A  and  m  denote  the  two  tuning  parameters:  A  controls  the  amount  of  smooth¬ 
ing  and  m  represents  the  number  of  continuous  derivatives  of  the  estimate.  Let  h  be  the  log  of 
the  spectral  density  and  h*  be  its  estimate  so  that  /„  =  exp(h«)  is  the  estimate  of  /.  Then,  by 
letting  A  be  a  function  of  the  sample  sise  n,  i.e.,  picking  a  sequence  A„  =  A(n),  which  decays 
to  0  slowly  enough,  certain  regularity  properties  of  hn  (and  hence  of  /„)  can  be  controlled.  In 
particular,  we  establish  that  if  A„  =  0{r»-r)  and  0  <  r  <  m/3,  then 

1.  supw  |hn(w)|  —  0,(1)  which  implies 

(a)  supw  /„(«)  =  0,(1)  and 

(b)  supw/-l(w)  =  0,(1) 

2.  sup  j/i|,(w)|3  =  0,(n1^*)  which,  together  with  condtion  (l.a),  implies  sup  |/i(«v)|*  = 

3.  E  /  |fc(u>)  —  /t»(w)|Jd<»>  =  o(l)  which,  together  with  condition  (l.a),  implies  /  |/(u>)  — 

/-MI’du.  =  0,(1). 

These  properties  can  be  seen  to  be  sufficient  for  the  conditions  described  at  the  end  of  the  previ¬ 
ous  section  to  hold,  and  thus  enable  the  results  of  this  paper  to  obtain  in  practical  applications 
where  the  parameter  of  interest  can  be  expressed  as  a  well-behaved  function  of  linear  functionals 
of  the  spectral  density. 
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1.3.4  Applications 

We  apply  this  resampling  methodology  to  the  proW»m  of  estimating  the  first  correlation  p  of 
the  process  by  the  sample  correlation  p  Since  p  and  p  can  be  expressed  as 

ratios  of  linear  functionals  of  the  spectral  density  and  periodogiam  respectively,  the  theoretical 
results  obtained  apply  in  this  case. 

The  simulations  studies  indicate  that  the  estimate  of  variability  obtained  by  cur  method  are 
roughly  as  good  as  the  asymptotic  approximations,  and  much  better  when  estimating  bias,  even 
though  the  asymptotics  were  calculated  assuming  complete  knowledge  of  the  model  generating 
the  data. 


1.4  Improved  estimators 

Suppose  we  observe  a  realization  of  size  n,  X£~l ,  of  a  Gaussian  stationary  process  with  spectral 
density  /,  and  we  estimate  a  general  functional  0/  =  0(f)  by  0Z. 

1.4.1  Circular  process 

Under  the  assumption  that  the  process  is  circular  the  periodogram  /*  is  a  sufficient  statistic 
and  we  can  consider  the  estimator  0'  =  E  which  has  smaller  risk  than  0Z  under  convex 

loss  by  the  Rao-Blackwell  Theorem. 

Under  th»  assumption  of  circularity  it  is  easy  to  calculate  0*  since,  conditional  on  the  pe¬ 
riodogram,  the  first  half  of  the  discrete  Fourier  transform  of  the  data  has  fixed  amplitude  and 
independent  uniform  U( 0, 1)  phases.  Thus,  by  sampling  iid  U( 0, 1)  phases  we  can  obtain  new 
samples  of  the  discrete  Fourier  transform  (DFT)  of  the  original  data  with  the  same  conditional 
(given  the  periodogram)  distribution.  Taking  inverse  DFT  we  can  obtain  new  samples  of  data 
Yq’*-1,  calculate  0,  for  each,  and  average  to  obtain  an  approximation  to  0*. 

1.4.2  Non-circular  process 

If  the  original  process  is  not  circular  the  periodogram  is  not  a  sufficient  statistic  but  we  may 
pretend  it  is  and  resample  as  in  the  circular  case  to  obtain  0*.  However,  this  mechanism  will  in 
general  introduce  bias,  i.e.,  EB\  ^  EBZ ,  but  the  amount  of  bias  is  often  neglibible  compared  to 
the  reduction  in  variance  so  that  0*  has  lower  mean  square  error  than  Bz. 

1.4.3  Applications 

The  motivating  idea  is  to  obtain  more  efficient  estimators  in  situations  where  efficient  estimators, 
e.g.,  MLE,  cannot  be  calculated  easily  but  there  exist  inefficient  estimators.  This  happens,  for 
example,  in  the  context  of  estimating  the  crossing  rate  of  level  a  of  a  Gaussian  stationary  process. 
In  this  case  the  MLE  (undet  the  Toeplits  model)  can  only  be  evaluated  in  the  case  a  =  0  in 
which  case  it  is  given  by  .5  —  w"1  arcsin(p),  where  p  is  the  first  sample  correlation  coefficient. 
On  the  other  hand,  the  empirical  crossing  rate  <j>°a  which  simply  counts  the  number  of  crossings 
and  divides  by  n  —  1,  can  be  calculated  for  all  a. 

The  simulations  indicate  that,  for  the  case  a  =  0,  the  MSE  of  E  is  roughly  equal 

to  that  of  the  MLE,  and  much  smaller  than  that  of  Furthermore,  it  is  worth  mentioning 
again  that  our  procedure  has  the  the  advantage  over  the  MLE  that  it  can  be  calculated  for  all 
ot  while  the  MLE  is  limited  to  the  case  a  =  0. 
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1.5  Algorithms 

Resampling  plans  aie  by  natuie  computet  intensive.  For  this  reason,  one  of  the  emphases  of 
this  research  was  to  implement  algorithms  which  could  be  driven  by  the  FFT.  To  generate  71- 
long  realizations  of  stationary  Gaussian  processes  with  covariance  matrix  En  we  adapted  an 
algorithm  to  simulate  circular  processes  of  length  2n  due  to  Jonas  (1984)  whose  n-long  segments 
have  covariance  E„. 

1.6  Related  research 

The  idea  of  using  the  data  to  specify  a  distribution  from  which  to  resample  was  advanced,  in 
the  iid  case,  by  Efron  (1979).  The  idea  of  resampling  Gaussian  stationary  process  with  spectral 
density  specified  by  the  data  to  study  sampling  bias  and  variability  was  set  forth  in  Ramos 
(1984).  Subsequently  Hurvich  and  Zeger  (1986)  and  Stine  (1986)  have  come  up  with  similar 
schemes  for  different  applications.  In  related  research,  methods  to  estimate  variability  in  the 
context  of  time-dependent  observations  have  been  suggested  by  Carlstein  (1986)  who  estimates 
it  on  the  basis  of  subseries  values  and  by  Freedman  and  Peters  (1984)  who,  in  the  original 
spirit  of  Efron  (1979),  suggest  fitting  a  long  autoregression  and  generating  new  psuedo-series  by 
resampling  the  residuals. 


1.7  Contents 

The  paper  is  organized  in  the  following  manner: 

1.7.1  Chapter  2 

In  this  chapter  the  background  and  numerous  auxiliary  results  are  presented.  Of  particular 
interest  are  results  concerning  the  magnitude  of  the  covariance  between  periodogram  ordinates 
and  their  logs,  Lemmas  2.16  and  2.17. 

1.7.2  Chapter  3 

In  this  chapter  we  establish  the  results  concerning  linear  functionals.  The  relevant  results  are 
Corollaries  3.5  and  3.9  which  establish  tha.  the  bootstrap  distribution  is  close,  in  probability, 
to  the  unknown  sampling  distribution  of  the  estimator. 

1.7.3  Chapter  4 

In  this  chapter  we  establish  the  properties  of  Wahba's  (1980)  procedure.  The  main  results  are 
Lemmas  4.3,  4.4,  4.8  and  4.7  which  establish  regularity  properties  and  convergence  in  Lj  norm 
to  0  of  the  integrated  mean  square  error. 

1.7.4  Chapter  5 

In  this  chapter  we  describe  the  Monte-Carlo  algorithm  for  simulating  realisations  from  Gaussian 
stationary  processes. 
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1.7.5  Chapter  8 

In  this  chapter  we  describe  applications  of  our  resampling  methodolgy  to  the  problem  of  esti¬ 
mating  the  sampling  bias  and  variability  of  the  first  sample  correlation  under  the  AR(1)  and 
MA(1)  models  and  result  of  simulations. 

1.7.6  Chapter  7 

In  this  chapter  we  examine  the  problem  of  finding  better  estimators  by  resampling.  An  applica¬ 
tion  to  estimating  the  crossing  rate  of  a  certain  level  is  discussed  and  results  of  simulations  are 
presented. 

1.7.7  Appendix  A 

All  proof  of  results  are  given  in  the  Appendix. 


Chapter  2 

Background 


In  this  chapter  we  state  the  basic  definitions,  assumptions,  and  properties  of  stationary  Gaussian 
processes  and  their  spectra;  we  also  define  the  discrete  Fourier  transform  (DFT)  and  the  peri- 
odogTam,  investigate  their  distributions  and  moments,  and  review  some  asymptotics.  Finally, 
we  briefly  discuss  circulants  with  two  examples  which  will  be  useful  later  on. 


2.1  Gaussian  stationary  processes 

In  this  section  we  introduce  Gaussian  stationary  process  and  investigate  properties  of  their 
covariance  structure  and  their  spectra.  We  also  describe  a  representation  of  such  processes  due 
to  Cramer. 

2.1.1  Covariance  structure 

In  what  follows  let  {Jf(}5^._QO  (or  more  simply  {Ai})  denote  a  real  and  discrete  time  series, 
that  is,  a  random  sequence  . . . ,  X_i,  Xo,  X\, . .  ■  .  In  this  paper  we  will  only  be  concerned  with 
time  series  which  are  stationary  and  Gaussian.  Intuitively  speaking,  stationarity  means  that  the 
distribution  of  the  X’s  does  not  change  with  time;  formally; 

Definition  2.1  {At}  is  a  Gaussian  stationary  sequence  if  all  its  finite  dimensional  distributions 
are  multivariate  r.  yrmal,  E(Xt )  =  n,  and  cov(AeA,+f)  —  c,,  for  allt,  s. 

Remark;  For  convenience  and  without  undue  loss  of  generality,  we  will  assume  throughout  this 
paper  that  n  =  0. 

Stationarity  implies  that  if  X$ ~l  =  (Xo,  ■  •  • ,  Xn-\ )'  is  a  realization  of  sise  n  of  a  stationary 
Gaussian  process  with  covariance  sequence  {c,},  its  covariance  matrix  En  is  Toeplitz  with  first 
row  c0, . . .  Recall  that  =  (<r;*)  is  Toeplits  whenever  ,  i.e.,  the  value 

of  any  entry  is  *  function  only  of  its  distance  to  the  main  diagonal.  Thus,  it  is  easy  to  see 
that  the  first  row  (or  column)  defines  any  Toeplits  E«  and  thus  it  is  convenient  to  write  En  = 
T(co, . . . ,  Cn_i)  where  T  maps  the  column  vector  (co, . . . ,  Cn_i)'  to  the  Toeplits  matrix  with  first 
row  (co, . . .  ,cn_t).  An  easy  consequence  of  definition  2.1  is  that,  for  all  n  and  s,  ~ 

1V(0,  £„)  where  E„  =  T(c0, . . .  ,c„_i). 

2.1.2  Spectrum 

A  dual  representation  of  the  covariance  structure  of  a  stationary  process  is  given  by  the  spectrum. 


7 


CHAPTER  2.  BACKGROUND 


8 


Definition  2.2  If  y)°°^  |c,|  <  oo  the  spectrum  of  the  process  can  be  defined  by 

OO 

/<*)=  E  c,  exp(i2irs  A), 

#  =  — OO 


(2.1) 


for  all  A. 

Assumption  2.1  The  covariances  of  all  processes  considered  in  this  paper  satisfy 

1.-00  |C.I  <  °°- 


The  characteristic  properties  of  the  spectrum  are: 

1.  /  is  clearly  periodic  with  period  1  so  we  can  limit  the  range  of  A  to  any  interval  of  unit 
length. 

2.  Since  {.Y,}  is  real-valued,  c,  =  c_,  and  hence  /(A)  =  /(  —  A)  which,  together  with  the  fact 
that  /  is  periodic,  implies  /  is  symmetric  about  1/2,  i.e.,  /( 1/2  -  x)  =  /(1/2  +  x). 

3.  /  and  {c,}  are  a  Fourier  transform  pair  and  hence 


/a-f  1 

/(A)exp(-t2i*A)dA. 


for  all  real  a  and  integer  s. 

4.  Brillinger  (1980,  p.  24)  shows  that  Definition  2.2  implies  f  >  0. 


(2.2) 


Lemma  2.1  It  f  >  0,  then  T(c0, ....  Cn-\)  >  0  for  all  n. 

Remark:  The  spectrum  determines  T(co, ....  cn_i )  for  all  n  so  that  we  will  write  Xg~l  ~  1V(0,  /) 

to  indicate  .YJ"1  ~  lY(0,T(co,...,cn.1)). 


Smoothness  of  the  spectrum 

It  is  easily  seen  from  (2.2)  that  a  constant  spectrum  implies  c,  =  0  for  «  ^  0,  i.e.,  the  { JYt }  are 
independent  and  the  process  is  memoryless.  On  the  other  hand,  if  the  spectrum  has  a  sharp 
peak  near  the  origin  the  c,’s  will  decay  slowly  and  we  can  think  of  the  process  as  having  a  long 
memory. 

In  the  next  few  paragraphs  we  briefly  discuss  in  what  way  “smoothness”  properties  of  the 
spectrum  and  conditions  on  the  covariance  sequence  {c, }  are  related  and  state  the  “smoothness” 
conditions  we  shall  impose  on  the  spectra  under  consideration  in  this  paper. 

“Smoothness”  of  a  function  /  is  usually  quantified  in  one  of  two  ways:  (i)  by  properties  of 
/,  namely  Lipschits  conditions  and  differentiability,  or  (ii)  properties  of  the  Fourier  coefficients 
{c,}  of  /,  namely  the  rate  of  decay  to  0  of  {c,}. 

Recall  the  definition  of  a  Lipschits  condition: 

Definition  2.3  /  €  Lip(a)  whenever  there  exists  a  smallest  constant  M/,a  that 

sup  \f(w  +  h)  -  f(u>)\  <  Mf,aea 


for  all  e  >  0. 
Remarks: 


1.  Mfa  is  called  the  a-modulus  of  continuity  of  /. 
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2.  We  will  write  Mj  for  Mj,\. 

It  is  clear  from  definition  2.2  that  ^  |c,  (  <  oo  implies  aupz  f(z)  <  co,  although  the  converse 
does  not  necessarily  hold.  In  general,  conditions  on  how  quickly  the  sequence  {c,}  converges  to 
zero,  differentiability  of  /,  and  Lipschitz  conditions  are  related  in  the  following  manner: 

Lemma  2.2  If  f(w)  =  ^^L_ao  c>  exp(i2rsw)  then 

<  oo  =>  sup|/<4)(w)|  <  oo  =.  /(k_l)  6  Lip(l) 

W 

for  k  >  1,  and  Assumption  2.1  implies  that  f  is  continous. 

Assumption  2.2  In  this  paper  we  shall  only  consider  processes  for  which  the  spectrum  f  6 
Lip(l)  and  f  >  0. 

2.1.3  Cramer  representation 

Throughout  this  paper  we  will  appeal  to  a  representation  of  Gaussian  stationary  processes  due 
to  Cramer.  Its  importance  lies  in  the  fact  that  it  expresses  the  values  of  the  process  as  a  simple 
function  of  the  spectrum  and  is  thus  useful  in  establishing  results  which  depend  on  properties 
of  the  spectrum. 

Theorem  2.3  (Cramer)  If  {Xt}  ~  N( 0,  /)  we  can  write 

Xt  =  /  exp(i2Ttfc>)F(ii/)dZ(u;)  (2.3) 

J-l/3 

where  F  is  any  integrable  complex  function  such  that  F( i*>)  =  F(-w)  and  |F(u>)|J  =  f{w),  and 
Z(w)  is  a  continuous  and  complex- valued  process  with  independent  Gaussian  increments  having 
the  following  properties: 

1.  EdZ{w)  —  0  for  all  w, 

2.  EdZ(w)dZ(X)  =  6u.xdu>, 

3.  EdZ(w)dZ(w)  =  0, 

4.  dZ{w )  =  dZ(-w). 

where  6Z  is  Kronecker's  6,  i.e.,  6t  =  1  if  x  =  0  and  0  otherwise. 

Remarks: 

1.  The  last  property  ensures  the  process  is  real-valued. 

2.  In  the  engineering  literature  F  is  called  a  transfer  function. 

3.  If  /  €  Lip(l)  it  will  be  convenient  to  take  the  transfer  function  F  €  Lip(  1 )  too.  In  other 
words,  let  F(w)  =  -J f(w)  exp(t£(u>))  where  6  6  Lip(l).  For  our  purposes  we  rarely  need 
to  consider  cases  where  0(w)  f-  0. 

2.2  The  main  construction  (MC) 

In  order  to  establish  the  main  results  of  this  paper  it  is  necessary  to  construct  two  processes 
{ }  ~  N( 0,  /)  and  {Y,}  ~  N( 0,  g)  defined  on  the  same  probability  space.  These  two  processes 
are  defined  in  such  a  way  that  they  are  to  be  “close”  (in  a  sense  to  be  made  precise  below)  to 
each  other,  subject  to  the  restriction  {Xj}  ~  N(0,  f)  and  {Yf}  ~  N(0,  g). 

First  we  state  the  assumption  on  the  transfer  functions: 

Assumption  2.3  Let  F ,  G  €  Lip(l)  be  measurable,  complex-valued,  non-vanishing  functions 
such  that  F(u>)  =  F(-w),  |F(ui)|l  =  /(u>),  G(w)  =  G(-w),  and  |C?(<*>) l3  =  s(u>). 
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The  main  construction:  Define  two  processes  {.Y,}  and  {Y,}  by 

fi/> 

Xt=  ej.p(i2rtw)F(ui)dZ(w)  (2.4) 

J- i/j 

and 

f 1/1 

Y,  =  /  exp(i2  rtu)G(u)dZ(u>)  (2.5) 

J-i/i 

for  all  t,  where  Z(u>)  is  as  in  Theorem  2.3,  and  F  and  G  satisfy  Assumption  2.3.  We  shall 
henceforth  refer  to  this  construction  as  the  Mam  Construction  (MC). 

Notice  that,  under  the  conditions  of  the  MC,  for  given  spectra  /  and  g  the  corresponding 
transfer  functions  F  and  G  are  defined  up  to  amplitude  but  not  phase  so  that  different  choices 
of  phase  functions  will  lead  to  processes  which  are  more  or  less  “close”  to  each  other.  If  we 
measure  “closeness”  by  E\Xt  -  Yt  |J,  it  is  then  straightforward  to  show  that 

Lemma  2.4  If  {Xt}  and  {Y,}  are  given  by  the  MC,  then 

E\Xt  -  Yt|3  =  J  |F(u>)  -  G(w)\1dui 

for  all  t,  and  is  minimized  when  the  phase  functions  of  F  and  G  are  the  same,  in  which  case 
E\Xt  -  Y,|*  =  / !  v/TT^T  —  s/g(wj\2du. 

2.3  The  DFT  and  the  periodogram 

In  this  section  we  define  the  discrete  Fourier  transform  (DFT)  and  the  periodogram  associated 
with  a  realisation  Xq~1.  We  also  examine  their  sampling  properties  when  Xq~1  ~  JV(0,/), 
/  €  Lip(l). 

2  3.1  Definition  of  the  DFT  and  periodogram 

Definition  2.4  Let  J*(w)  =  n-1^’  Xte~'lTwt  denote  the  discrete  Fourier  transform  of 

yn  —  1 

AQ 

Definition  2.5  Let  /*(<*))  =  |7*(«)|J  be  the  associated  periodogram. 

Remits: 

1.  J*(w),  /*(w)  are  periodic  with  period  1. 

2.  If  Xq~1  is  real-valued  it  follows  that  J*(w)  =  /*(1  -u>).  Hence  /*(w)  =  /*(1  -  u<)  and 

We  will  write  /*  without  the  argument  u>  to  denote  the  column  vector  with  entries  J*(k/n), 
k  =  0, . . .  ,n-  1,  and  similarly  /*  is  the  vector  whose  components  are  |/*(lk/n)|J,  k  =  0, . . . ,  n- 1. 
These  are  just  the  DFT  and  the  periodogram  evaluated  at  the  set  of  Fourier  frequencies  Jfc/ n, 
k  =  0, . . . ,  n  —  1.  We  will  refer  to  these  vectors  simply  as  the  DFT  and  the  periodogram. 
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2.3.2  Convolution  and  kernels 

Since  many  of  the  results  we  derive  below  are  expressed  in  terms  of  convolutions  of  a  function 
/  and  a  kernel,  either  Fejer’s  or  Dirichlet’s,  we  state  the  the  relevant  definitions  and  properties 
below.  In  this  section  /  and  g  denote  arbitrary  periodic  functions  (not  only  spectral  densities). 

Definition  2.8  If  f  and  g  art  periodic  define  the  convolution  f  *  g  of  f  and  g  by 


f*g(z)  =  J  f{*-t)g( t)dt  =  J  f{t)g(x  -  t)dt 


Remarks: 


1.  We  will  usually  take  a  =  0  or  a  =  -1/2. 

2.  Notice  that  if  either  /  of  g  are  differentiable,  their  convolution  is  differentiable  and 

(/**(*))'=  f  f'(z-i)g(t)dt=  f  f(t)g'(z-t)di. 

Jo  Jo 


In  this  paper  we  will  use  the  following  notation:  let  Lp  denote  the  space  of  periodic  (almost 
everywhere)  functions  which  are  Lebesgue  integrable  to  the  p,h  power  with  norm 

i.e.,  /  €  Lp  if,  and  only  if,  ||/ ||p  <  oo.  If  /  and  g  are  in  Lp  we  will  occasionally  write  dp(f,  g)  = 

11/  ~  9\\p- 

Lemma  2.5  Let  fyg  g  L\,  g  >  0,  and  ||p||i  =  1.  Then 

inf /(w)  <  inf  /  •  g(a>)  <  sup  /  *  g(u>)  <  sup/(u»). 

^  W  u  y 


In  the  following  lemma  5  denotes  either  (i),  the  set  of  periodic  continuous  functions  with 
sup  norm,  or  (ii),  Lp. 

Lemma  2.8  (Butzer  and  Nestel,  1971).  Let  f  S  5  and  g  £  L\.  Then  f  •  g  £  S  and  \\f  *  g\\s  < 

H/iisIMli. 

Remark:  If  ||yj|i  =  1  we  can  apply  Lemma  2.6  to  the  derivatives  of  /  to  show  that  the  result  of 
convolving  /  and  g  is  smoother  than  /. 

We  now  state  definitions  and  properties  of  the  Dirichlet  and  the  Fejer  kernels. 


Definition  2.7  For  n  >  Q  let 


denote  Dirichlet 's  kernel  and  let 


DnW  = 


sin(nxA) 

sin(irA) 


*n(A)  = 


sin’fnwA) 

nsin3(vA) 


denote  Fejer '»  Kernel. 
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These  kernels  arise  naturally  when  approximating  /  by  finite  “sums”  of  its  Fourier  coefficients 
{c*}.  Here  the  word  “sum”  may  refer  to  different  types  of  sums,  e.g.  Cesaro.  Each  of  these 
“sums”  can  be  expressed  as  the  convolution  of  /  with  a  particular  kernel. 

The  following  are  standard  results  from  approximation  theory  (cf. ,  e.g.,  Butser  and  Nessel, 
1971). 

Lemma  2.7  (Butzer  and  Nessel,  1971)  If  D„  and  Kn  are  the  Dirichlet  and  Fejer  kernels  re¬ 
spectively,  then 

n 

T.  exp{t2rfcA)  =  D}n  +  i(A),  (2.6) 

k  =  -  n 

n  - 1 

exp(  -ix(n  -  1)A)  ^  eip(i2xk\)  =  Dn(\)  (2.7) 

h=0 

and 

y  (l  -  l—^j  exp(i2**A)  =  Kn(X).  (2.8) 

Remark:  Since  Kn  >  0,  Lemma  2.7  implies  ||ffn||i  =  1. 

Lemma  2.8  If  f  €  Li  and  c*  =  fQl  f( A)  exp(- i2ts A)  <fA,  then 

n 

y  ckcxp{i2xkx)  =  Dj„+i  *  f(i) 

k~-n 

and 

^2  Ck  (l  “  =  Kn  *  f(x). 

k=-n  '  U  ' 

Lemma  2.9  (Butzer  and  Nessel,  1971)  If  f  €  Lip(l),  then 

sup  I  Jr*  *  /(w)  -  /(<*>)!  =  Af/0(logn/n). 


uniformly  in  f. 

The  next  is  an  auxiliary  lemma  we  will  use  often. 
Lemma  2.10  If  Dn  is  the  Dirichlet  kernel,  then 


/«-f  l 

Dn(x  -  t)Dn(x  -  s)dx  =  Dn(s 


-<)■ 


for  all  a.  If  in  addition  f  €  Lip{  1)  is  symmetric  and  periodic,  then 

f  •+ 1 


/•-M 

{/(*)  -  f(')}Dn(x  -  t)Dn(x  -  s)dx  =  MfO(\ ogn). 


uniformly  in  /  and  hence 
*+i 


i: 


f(*)Dn(x 


-  t)Dn(x  -  s )dx  -  Dn(s  -  t)f(s)  +  MfO(logn). 


uniformly  in  f  for  all  s,t,  and  a. 
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2.3.3  Distribution  and  moments  of  the  DFT 

Let  F„  denote  the  matrix  with  entries  n-1/5  exp(-i2irjfc/n),  j,  k  =  0,  ...,n  -  1.  Fn  is  called  the 
Fourier  transform  matrix  since  it  projects  a  vector  onto  its  DFT,  e.g.,  J*  =  FnXZ~l.  Notice 
that  Fn  is  symmetric  and  unitary,  i.e.  Fn  1  =  Fn.  Fn  should  not  be  confused  with  F,  the 
transfer  function. 

If  we  write  Fn  =  Cn  -  iSn  for  real-valued  matrices  Cn  and  Sn,  we  can  establish  the  following 
result  concerning  the  distribution  of  J*: 

Lemma  2.11  If  =  F„X£~l  and  X^1  ~  N(0,En),  then 


CnEnCn 

~F, v  £nC„ 


-CnEnSn  \\ 

SnSnSn  j)' 


and  F|7*(u/)|J  =  EI*(w)  =  Kn  .  f(u>)  for  alia). 

The  following  lemma  gives  the  variances  and  covariances  between  real  and  imaginary  parts 
of  the  DFT. 


Lemma  2.12  If  X£ ~l  ~  N( 0,  /)  for  f  6  Lip(l)  and  0  <  I,  m  <  [n/2],  then 

2 FK{7*(//n)}S{J*(m/n)}  =  5m_r(  1  +  6t  +  6n/1_,)Kn  »  f{m/n)  +  MfO( logn/n), 
2F3{^(l/n)}3{7„’(m/n)}  =  1  -6,-  6n/i_,)Kn  .  f(m/n)  +  MfO(\ ogn/n), 

and 

ES{j;(l/n)}S{J:(m/n)>  =  (1  -  6m)(  1  -  6 - n)M,0{\ ogn/n) 

uniformly  in  f  £  Lip{  1). 


2.3.4  Distribution  and  moments  of  the  periodogram 

Distribution  of  periodogram  ordinates 

We  saw  above  (Lemmas  2.11  and  2.12)  that  the  real  and  imaginary  parts  of  J*(k/n)  have 
bivariate  normal  distribution  with  variance  approximately  proportional  to  Kn  »  f(k/n)/2  and 
correlation  p  —  MfO(log  n/n).  Lemma  2.11  implies  that  twice  the  periodogTam  ordinate  I*(k/n) 
is  distributed  like  WKn  *  /(k/n)  where  W  =  X 3  +  (pX  +  \/l  -  and  X  and  Y  are  iid 

N( 0, 1)  variates.  The  following  lemma  establishes  an  alternate  representation  for  W  which  will 
make  manipulations  easier. 

Lemma  2.18  If  W  =  X*  +  (pX  +  y/l  -  p3Y)i  and  X  and  Y  are  ixd  1V(0, 1),  then  W  is 
distributed  as 


w  =  (\  +  p)ui  +  [\ -  p)ul 

=  (£f*  +  ^j)(l  +p€) 


for  Ui,Uj  iid  N(0, 1).  The  two  factors  on  the  r.h.s.  are  independent  and  \0\  <  1. 
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Covariance  between  periodogram  ordinates 

Below  we  establish  the  main  result  of  this  section  which  calculates  the  covariance  between 
periodogram  ordinates  of  realisations  X£~l  and  V0n_l: 

Lemma  2.14  If  Xq~1  and  V0n_1  are  given  by  the  MC  with  transfer  functions  F  and  G,  then 
cov(Fn(\ ),/»(#*))  =  n-^J  F(u)G(Z)Dn(u,-\)Dn(u>-»)dw 

+  n-3  |  J  F(u/)G(ui)Dn(w  —  A)Dn(w  +  fi)dw 


In  the  case  F  —  G,  i.e.,  A',  =  K,,  lemma  2.14  leads  t"  the  following  known  result  (cf.,  e  g., 
Hassminski  and  Ibragimov,  1986): 

Corollary  2.15  // AJ'1  ~  JV(0,  /) ,  then 


cov(i;(A  ),/;(/*)) 


=  n~3  J  J  f(u>)Dn(w  -  \)Dn(u>  -  p)<b 

+  n"3  j  f  f(ui)Dn(w  -  A )Dn(u>  +  n)du 


(2.9) 


It  is  well  known  that  the  components  of  /*(fc/n),  k  =  0 . [n/2]  become  increasingly 

uncorrelated  as  n  — •  oo.  Since  spectral  estimates  usually  depend  only  on  these  components,  the 
properties  of  these  estimates  may  depend  on  the  rate  at  which  the  correlation  between  these 
components  decays  to  0.  For  this  reason,  it  is  useful  to  evaluate  (2.9)  explicitly  when  A  and  n 
are  different  Fourier  frequencies.  We  then  have  the  following 

Lemma  2.16  // AJ~l  ~  JV(0,/),  /  €  Lip(l),  and  0  <  I  <  m  <  [n/2],  then 


cov(J2  (*/»).£  ("*/*>))  =  M;sO(logn/n)3. 

uniformly  in  f. 

We  conclude  this  section  with  a  result  on  the  magnitude  of  the  covariance  between  the 
logs  of  periodogram  ordinates.  This  result  will  prove  useful  when  establishing  properties  of 
the  particular  spectral  estimate  (Wahba,  1980)  used  in  our  simulations  since  this  procedure 
essentially  amounts  to  taking  a  weighted  average  of  the  log  of  periodogram  ordinates. 

Lemma  2.17  //A"-1  satisfies  the  conditions  of  Definition  2.3  and  0  <  1  <  m  <  [n/2],  then 

cov(log  /*(I/n),  log  7*(m/n))  =  Af30(log  n/n)3.  (2.10) 


uniformly  m  /. 


2.4  Circulants 

In  this  section  we  introduce  circulants,  describe  some  of  their  properties,  and  discuss  in  some 
detail  two  examples  of  circulants  which  will  appear  later  on. 

Definition  2.8  A  matrix  E'  =  (<r*y),  k,j  =  0,  ...,n  -  1,  is  a  circulant  of  dimension  n  if 
<D,j  =  <re(ifc  —  j)  for  some  function  ffc  with  period  n. 
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Recall  (cf.,  e.g.,  Brillinger,  1979)  that  the  eigenvalues  of  E‘  are  given  by 

n-l 

=  ^2<rc(k)eip(-t2iek}/n), 
k=0 

j  =  0, n  —  1,  and  the  corresponding  right-hand  eigenvectors  are  simply  the  rows  (or  columns) 
of  the  Fourier  transform  matrix  Fn.  In  fact,  it  is  well  known  (cf.  e.g.  Brillinger,  1979)  that  E' 
is  a  circulant  if,  and  only  if,  it  can  be  written  as 

Zcn  =  FnA'nTr  (2.11) 

where  A'  is  diagonal  with  entries  AJ.  The  importance  of  (2.11)  is  that  the  eigenvectors,  the 
columns  of  Fn,  are  the  same  for  all  circulants  of  dimension  n.  Notice  that  whenever  E'  and  A' 
are  both  real,  E*  is  necessarily  symmetric  (and  hence  Toeplits)  and  E‘  =  FnA^Fn. 

It  is  easy  to  see  that,  given  the  eigenvalues,  we  can  obtain  the  entries  of  EJ;  by  taking  the 
inverse  DFT 

n-l 

»'(*)  =  ~  H  Ar  exP (»2»>*/n).  (2.12) 

n 

1=0 

2.4.1  Example  1:  A  circulant  approximation  to  En 

Let  Xq  -1  be  a  realisation  of  sise  n  from  JV(0,  /),  /  >  0,  and  let  E«  denote  its  covariance  matrix, 
i.e.,  Xg~l  ~  !V(0,£n),  E„  =  T(cq,  . . .  ,cn_1).  If  we  let  J*  =  FnX  denote  its  DFT,  it  follows 
from  Lemma  2.11  that 

EJ^lf  =  FnZnK  =  A„. 

so  that 

=  F»A,Fn  =  FnKFn. 

where  A*  is  not  necessarily  diagonal  or  real.  However, 

Lemma  2.18  The  diagonal  entries  of  An  ore  given  by 

n-l 

Ay  =  c‘  0  -  i*l/n)  exp(-«2)rt>/n)  (2.13) 

t=-n+l 

n-l 

=  (M1  -  l‘l/»») +  c„_,|t|/»)exp(-«2irtj7»)  (2.14) 

«=o 

Remark:  The  Ay  =  Kn*  f{j /n)  are  the  Cesaro  means  of  /  evaluated  at  the  Fourier  frequencies 

>/»,  j=  1. 

In  order  to  construct  a  circulant  approximation  to  En,  define  A£  =  diag(An)  to  be  the 
diagonal  matrix  whose  entries  coincide  with  the  entries  in  the  main  diagonal  of  A*  and  let 

Z'n  =  Fn\'nK  =  F;KFn. 

Lemma  2.12  shows  that  the  off-diagonal  entries  of  A„  are  0(logn/n)  so  that  £'  is  a  circulant 
approximation  to  En. 

By  construction  it  is  clear  that  E'  is  a  circulant  since  it  is  in  the  form  of  (2.11).  The  entries 
of  A*  are  given  by  (2.14)  and  hence  the  inversion  formula  (2.12)  gives 

cj  =  c4  (1  -  fc/n)  +  cn_»k/n.  (2.15) 

for  k  =  0, . . . ,  n  -  1.  Notice  that  the  approximation  is  better  for  entries  closer  to  the  main 
diagonal. 
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2.4.2  Example  2:  “Circularizing”  processes 

Let  f(u/)  =  ^2,  c,  exp(»2rws)  be  a  spectrum  and  consider  the  following  2n  x  2n  circulant 

s3»  =  T(c0,...,cn_,,0,cn_i,...,ci). 

Lemma  2.19  //  £}„  -  FinN \nF j„,  then  Ajn  is  diagonal  with  entries 

n  —  1 

AJ  =  ^2  C,exp(i2irt;/2n) 

t  =  —  n  + 1 


for  j  =  0, . . . ,  2n  -  1. 

Remarks: 

1.  The  eigenvalues  A'  =  Djn+ 1  *  f(j/2n)  are  the  partial  sum  of  the  Fourier  series  of  / 
evaluated  at  the  Fourier  frequencies  ji/2n,  for  j  =  0, . . . ,  2n  —  1.  Since  the  Dirichlet  kernel 
is  not  positive  everywhere  it  may  be  that  some  AJ’s  are  negative  even  when  /  >  0,  and 
hence  that  £|n  is  not  necessarily  positive  definite. 

2.  If  Aj  >0  for  all  j,  then  £|B  is  positive  definite  and  hence  is  the  covariance  matrix  of 
a  circular  process  with  period  2n.  This  process  has  the  property  that  any  n  adjacent 
components  of  a  realisation  are  distributed  jV(0,£„)  where  E„  T(co, . . . ,  cn_i ). 

The  last  property  will  be  useful  for  fast  simulation  of  realisations  from  7V(0,  £„)  since  the  FFT 
can  be  used  to  generate  realisations  of  processes  which  have  circulant  covariance  matrices.  For 
this  reason,  it  will  be  important  that  for  large  enough  n,  £jn  be  positive  definite,  i.e.,  that  it  be 
a  covariance  matrix.  This  is  guaranteed  by  the  following 

Lemma  2.20  //£|c*|  <  oo  and  f  >  0,  then  there  exists  N  such  that  is  positive  definite 
for  n  >  N . 


Chapter  3 


Linear  functionals  of  the 
spectrum 

3.1  Introduction 

Recall  the  motivating  problem  stated  in  the  introduction:  upon  observing  a  realisation  Xq  ~ 
N(0,f) ,  /  unknown,  procedure  0,  is  used  to  estimate  a  functional  9f  of  f  and  it  is  of  interest 
to  assess  the  sampling  distribution  of  0,. 

One  possible  solution  to  this  problem  is  the  following:  use  some  procedure  /„  to  estimate  / 
and,  by  simulating  many  realisations  from  a  process  {F,}  whose  conditional  distribution  given 
/„  is  N(Q,  /„),  determine  the  (conditional)  sampling  distribution  of  0,  and  use  it  to  estimate 
that  of  0,.  An  obvious  question  is:  under  what  conditions  on  0,  0,  /,  and  /„,  can  the  sampling 
distribution  of  9Z  be  consistently  estimated  by  the  (conditional)  sampling  distribution  of  0y? 

In  this  chapter  we  investigate  the  answer  to  this  question  in  the  case  where  the  parameter  of 
interest  is  a  Linear  functional  of  the  spectrum  and  its  estimate  is  the  corresponding  functional 
of  the  periodogram.  The  main  result  in  this  chapter,  Corollary  3-5,  establishes  that  in  this  case, 
under  certain  conditions  on  /„  and  0/,  C(9Z)  and  £(0y|A’o-1)  are  close  with  high  probability. 
This  result  is  then  extended  to  well-behaved  functions  of  linear  functionals  in  Corollary  3.9.  For 
example,  we  shall  be  interested  in  assessing  the  sampling  distribution  of  estimators  of  normalised 
linear  functionals  such  as  the  ratio  of  two  linear  functionals. 

Our  results  will,  of  course,  depend  on  properties  of  the  particular  estimator  of  the  spectral 
density  used.  In  Chapter  4  we  show  that  the  estimator  of  the  spectral  density  (Wahba,  1980) 
which  we  use  in  our  simulations,  will  have  the  required  properties  whenever  the  smoothing 
parameter  is  properly  chosen. 

Some  of  the  theory  required  for  establishing  these  results  can  be  found  in  §A.2.1. 


3.2  Sampling  distribution  of  estimators  of  linear  func¬ 
tionals 

3.2.1  Efficient  estimators 

Assume  we  observe  a  realisation  Xq  ~l  ~  1V(0,  /),  /  unknown,  and  we  estimate  the  parameter 

0/  =  f  0(w)/(i^)dto>  (3.1) 
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by 


9*  =  I  «(-)£ 


(3.2) 


for  a  teal  and  measurable  function  9.  There  is  no  loss  of  generality  in  assuming  that  9  is  even 
and  periodic  and  we  shall  hereafter  do  so. 

Linear  functionals  include  a  wide  class  of  pax  eters,  e.g.,  covariances  and  spectral  moments. 
In  addition,  they  have  nice  properties:  Hasminskii  and  Ibragimov  (1986)  establish  (under  weaker 
conditions  on  /  and  9  than  those  we  assume  in  this  paper)  that  9Z  is  asymptotically  efficient, 
and 

£  (>/S(«,  -9,))  -  C(U)  (3.3) 


where  U  ~  N  (0,  2  f  0(w)3/(u>)Jdw) . 


3.2.2  Resampling  strategy 

We  saw  above  that  \/n(6z  -  6/)  has  a  Limiting  distribution.  For  this  reason  it  will  be  easier  to 
deal  with  the  sampling  distribution  of 

Un  =  VZ(9z  (3.4) 


rather  than  that  of  9Z . 

The  idea  behind  our  resampling  methodology  is  to  estimate  /  by  /„  and  estimate  the  distri¬ 
bution  of  L7„  by  the  conditional  distribution  given  /„  =  g  of 

V»  =  VZ(9,  -  9,),  (3.5) 


where  \X^~: l)  =  N(0,  /„). 

3.2.3  Results 

The  results  of  this  section  are  obtained  by  showing  that  there  exist  random  variables  Un  and 
V",,  which  are  doae  with  high  probability  and  such  that  £(Cfn|A'J-1)  =  £(17n)  —  £(U„)  and 
£(V^|XJ-1)  =  £(V^|JfJ_l),  i.e.,  Un  -  Vn  =  op(l).  Corollary  3.5  then  establishes 

1.  dP  (£(ff,|.Xy-1)f£<?.|Xy-1))  =o,(l)  tnd 

2.  dP(c(U),C(Vn\XZ~1))  =or(l). 

where  dP  denotes  the  Prohorov  distance  (cf.,  e.g.,  JA.2.1).  In  other  words,  the  distribution  of 
Un  is  approximated  by  the  “bootstrap”  distribution. 

In  the  next  section  we  construct  these  random  variables  and  show  that  they  indeed  have 
the  desired  properties.  To  do  so,  we  proceed  in  stages:  first,  we  assume  /„  =  g  is  given  and 
deterministic  and  use  the  MC  to  construct  processes  {.Xi}  ~  W(0 , /)  and  (T)}  ~  AT(0,  j)  and 
define  Un  and  Vn  as  in  (3.4)  and  (3.5);  second,  we  bound  the  distance  between  the  distributions 
of  Un  and  Vn  in  terms  of  properties  of  /  and  g;  third,  we  extend  the  results  to  account  for  the 
fact  that  /„  is  really  a  random  variable  by  determining  regularity  conditions  on  /„  so  that  the 
required  properties  of  deterministic  g  hold  for  random  /„  with  high  probability. 
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3.3  Distances  between  distributions  of  estimators 

3.3.1  Deterministic  case 

If  /«  —  g  »*  deterministic,  the  basic  idea  is  to  study  how  “close”  the  distributions  of  Un  and 
Vn  are  as  a  function  of  how  “close”  the  respective  spectra  /  and  g  are.  It  turns  out  (Theorem 
3.3)  that  if  the  two  processes  {X(}  —  N(0,f)  and  {y,}  ~  N(Q,  g)  are  given  by  the  MC  and  we 
measure  the  distance  between  the  distributions  of  Un  and  Vn  by  E\Un  —  Vn|,  this  distance  is 
essentially  bounded  by  the  Lj  distance  between  /  and  g. 

Before  we  go  on  we  need  to  establish  two  auxiliary  results.  The  first  one  concerns  the  bias 
of*,. 

lemma  3.1  If  9  is  bounded,  9X  is  given  by  (3.2),  f  6  Lip(l),  and  X*~l  ~  N(0,f),  then 

E{9Z  -  9/)  —  M)  Oflog  n/n) 


uniformly  m  f . 

The  second  one  concerns  the  covariance  between  9X  and  <£,  =  f  <j>(ui)I*(w)du/,  an  estimate 
of  <f>}  =  f  <fr(u>)g(ui)du> . 

Lemma  3.2  //  X£~l  ~  N(0,f)  and  yon_l  ~  W(0,p)  are  given  by  the  MC  with  transfer  func¬ 
tions  F  and  G,  and  9,d>  €  Lip(l)  are  eoen  and  periodic,  then 

nco v(0„,<£,)  =  2  J  9(ui)d>(u/)f(u>)g(u)cL >  +  MF^  sup  |FG|0(log  n/n)  +  sup  |/g|0(logJ  n/n) 

uniformly  in  F  and  G. 

We  are  now  able  to  establish  the  establish  the  main  result  of  this  section,  Theorem  3.3,  which 
bounds  the  distance  between  the  distribution  of  Un  and  Vn  by  E\Un  -  Vn\  where  Un  and  Vn  are 
induced  by  processes  {X,}  and  {Ti}  given  by  the  MC. 

Theorem  3.3  //  Xj_l  ~  N(Q,f)  and  yon_l  ~  IV(0,g)  ore  given  by  the  MC  with  transfer 
functions  F  and  G,  9  €  Lip(l)  is  even  and  periodic,  and  Un  and  Vn  are  given  by  (3-4)  and 
(3.5),  then 

E'\Un-Vn\  <  E]Un  -V„|J 

=  2je(w),]f(w)-g(w)\idu,  (3.6) 

+  Mi?3Sup|FG|0(logn/n)  +  (sup  \fg\  +  M)  +  Af’)0(log3  n/n) 

uniformly  in  F  and  G. 

This  result  will  enable  us  to  show  that  Un  -  Vn  =  0,(1)  when  g  is  no  longer  deterministic  but 
a  consistent  (in  integrated  mean  square  error)  estimate  /„  with  certain  regularity  properties. 
Consistency  ensures  the  integral  in  (3.6)  to  be  0,(1)  while  the  regularity  conditions  ensure  the 
remaining  terms  in  (3.6)  to  be  o,(l). 

3.3.2  Random  case 

Extending  the  results  of  the  previous  section  to  the  random  case  where  }n  is  random  will  follow 
upon  showing  that  the  estimate  of  the  spectral  density  based  on  the  original  sample  XJ-1  and 
the  corresponding  estimate  of  the  transfer  function  satisfy  certain  properties  which  we  will  now 
explore. 
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Let  /»  denote  the  estimate  of  the  spectral  density  which  the  observed  data  X ™  1  induce. 
We  will  use  the  MC  with  F  =  G  =  and  Z(w)  independent  of  Xq~1  to  generate 

realizations  of  two  processes  {Jtt}  and  {y,}  whose  conditional  distributions  given  Xr1  (or 
equivalently  /„)  are  N(0,  f)  and  N(0,fn).  These  processes  in  turn  induce  random  variables  t/„ 
and  Vn  with  the  property  that  C(Un\X^~l)  =  £(</„)  =  C(Un)  and  £(V„|A-J-1)  =  £(Vnl*r1)- 
We  will  show  that  under  regularity  conditions  on  fn  these  two  distributions  are  close  with  high 
probability. 

Theorem  3.4  Let  /„  be  a  spectral  estimate  which  depends  on  Xq  -l  and  satisfies  the  following 
properties: 

!■  11/ -/„||i  =<*(!) 

2.  P(infw  /„(*)  >  0)  -  1. 

3-  supw  fn(w)  =  0,(1) 

4-  *up„  /„_1(w)  =  Op(l) 

5 ■  *«P«  /AM1  =  of.(n/login). 

If  0  €  Lip(l)  is  even  and  periodic,  then 


Un-vn  =  0,(1). 

Corollary  3.5  Let  U  ~  Af(0,2  f  0(a))3  f(u>)3du) .  If  the  conditions  of  Theorem  3-4  hold,  then 

1.  dP  (tlUJ.TfKIXf-1))  =  .y(l)  and 

2.  dp  (C(U),  £(Vn|X"-1))  =o,(l). 

Remark:  Item  (1)  above  shows  that  the  sampling  distribution  of  U„  is  consistently  estimated 
by  the  “bootstrap”  distribution. 


3.4  Functions  of  linear  functionals 

In  this  section  we  extend  the  results  of  the  previous  section  to  well-behaved  functions  of  linear 
functionals  of  the  spectrum,  i.e.,  the  case  where 

is  used  to  estimate 

Wn . «/,)  =  M»/) 

for  some  function  A  which  is  well-behaved. 

3.4.1  Limiting  distribution 

In  this  section  we  establish  that  the  limiting  distribution  of  well-behaved  functions  of  linear 
functionals  is  multivariate  Gaussian.  First,  we  show  that  the  limiting  distribution  of  vectors  of 
linear  functionals  is  multivariate  Gaussian. 
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Lemma  3.8  Let  f ,  6}  6  Lip(l),  j  =  1, . . . ,  p,  and  lei 

be  a  vector  of  estimates  of 

*/=(«/! . «/,)' 

where 

9t,  =  j  6, (u,)rn(u,)d(u>) 

and 

9fj  =  y 

/or  Then,  the  limiting  distribution  of  y/n(8z  —8f)  is  A'(0,E)  where 

The  following  theorem  can  then  be  used  to  establish  asymptotic  normality  of  estimators 
well-behaved  functions  of  linear  functionals. 

Theorem  3.7  (Rao,  1972)  IfTn  =  (Tln, . . . ,  Tk„)  is  a  k-dimensional  random  variable  such  that 
the  asymptottc  dis.ribution  of 

/  ^n(Tln  -  6X)  \ 

V  / 

tj  k-variate  normal  with  mean  zero  and  covariance  matrix  E  =  (<rl;)  and  h  is  a  function  of  k 
variables  which  is  totally  differentiable,  then  the  asymptotic  distribution  of 

Mh(Tin,-.-,Tkn)-h(6l...8k)) 


is  normal  with  mean  zero  and  variance 


”(#)  = 

«  j 


dh  dh 

WtW, 


provided  v(0)  >  0,  6  =  (0j, . . . ,  9k)' . 

3.4.2  Results 

In  this  section  we  extend  the  result  of  Theorem  3.4  to  well-behaved  functions  of  linear  functionals 
of  the  spectrum. 

Theorem  3.8  Let  /„  be  an  estimate  of  the  spectrum  f  which  depends  on  X£~l  and  satisfies 
the  assumptions  of  Theorem  S.f.  Let  X*~l  and  Y"~l  be  realizations  of  two  processes  given  by 
the  MC  with  transfer  functions  F  =  ( / )l/,J  and  G  =  (/n)1^1,  and  let  8}  £  Lip(l),  j  =  l,...,p 
be  even  and  periodic.  If 

K  =  y/n(h(0t)  -  h(e,)),  v;  =  MW,)  -  WfJ), 

and  h  has  continous  derivatives  in  some  neighborhood  Ngf  of  8j,  then 

u;  -  v;  =  0,(1). 
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Corollary  3.9  Assume  the  conditions  of  Theorem  3.8  hold,  and  let  U"  be  a  random  variable 
whose  distribution  is  the  same  as  the  limiting  distribution  of  U° .  Then 


1.  aP  (C(U-),C(K I.Y"-1))  =  o„(l)  and 

2.  dp  (C(U-),C(V;\XZ-1))  =op(i). 


3.4.3  Application:  normalized  functionals 


We  illustrate  the  results  of  the  previous  section  with  an  example  which  we  use  later  on  in  our 
simulations. 

Theorem  3.7  can  be  used  to  establish  asymptotic  normality  of  estimates  of  normalized  linear 
functionals,  i.e.,  estimates 

f0(u,)mv)du,  ez 

*  JiZHdu,  ||/«||t-  1  ' 


of  functionals 

/ 6(u>)f(w)<bi  _  e, 
ff(»)du,  ll/lli  * 

The  class  of  normalized  functionals  includes  the  correlations 


(3.8) 


f  cos(2ru/)f(u/)dw 
f  f(u)du 


and  the  normalized  spectral  moments. 

A  corollary  to  Theorem  3.7  gives  the  form  of  the  limiting  distribution  of  normalized  linear 
functionals. 

Corollary  3.10  If  Un  =  \/n(8l  —  8)),  then  its  limiting  distribution  ts 


N  0 


,2(8))'  j (9(u>)/8)  -  \)' }'(*)' <L, 


where  /*  =  //li/Hi- 

Corollary  3.9  implies  that  the  sampling  distribution  of  Un  can  be  consistently  estimated 
by  our  resampling  methodology.  In  Chapter  6  we  give  the  results  from  simulations  where  the 
parameter  of  interest  is  a  normalized  functional,  the  first  correlation. 


Chapter  4 

Estimating  the  spectral  density 


4.1  Wahba’s  estimate 

In  this  chapter  we  introduce  Wahba’s  (1980)  estimate  and  establish  that  it  has  the  properties 
which  will  allow  us  to  apply  the  results  of  Chapter  3  when  it  is  used  as  the  spectral  density  to 
generate  the  bootstrap  samples. 

The  procedure  estimates  the  log-spectrum  of  a  stationary  process  with  positive  and  differ¬ 
entiable  spectrum  /  €  Lip(l).  In  this  case  the  log-spectrum  h  6  Lip(l)  and  has  a  Fourier  series 
expansion,  i.e., 

no 

h(u)  =  log/(w)  =  E  h»  exp(i2*wfc). 

kzz-oo 

We  will  assume,  in  addition,  that  |6»|  <  oo. 

Now,  suppose  we  observe  a  realisation  of  length  2n  from  a  Gaussian  stationary  process  with 
spectrum  /,  i.e.,  we  observe  Jfgn~l  ~  JV(0, /).  The  proposed  estimate  of  h(w)  is 

mw)=  e  ■llj;;)i.«p(,2^).  (4i) 

v=-n  '  ' 


where 

1  - 

6vn  =  2 -  E  YknCTp(-i2iruk/2n),  (4.2) 

k=~n+l 

Ykn  =  log/Jw(*/2n)  +  7,  and  7  =  .57721  ...  is  Euler’s  constant. 

From  (4.1)  and  (4.2)  it  is  easy  to  show  that  hn  is  a  smoothed  version  of  the  log-periodogram 
since  it  can  be  rewritten  as 

n 

M")  =  E  *»«W»,m,B(w  -  */2n)  (4.3) 


where 


exp  (i2rvw) 

1  4-  \(2rv)3m 


is  real,  periodic,  not  necessarily  positive,  and  such  that 


(4.4) 


n 

(w )dw  =  E  WVm.n(w  -  fe/2n)  =  1. 

k=  -  n  + 1 
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4.2  Properties  of  the  log-periodogram 

In  order  to  study  the  properties  of  (4.3)  it  is  convenient  to  first  assess  the  properties  of  the 
triangular  array  Y*n,  k  =  — n  +  1, . . . ,  n.  We  can  write 

=  log  (ffjn  *  /(*/2n)Vjn)  +  y 
=  \og(K3n*f(k/ 2n))+f*„ 
where  f*n  =  log  Vkn  +  y,  and  establish  the  following 

Lemma  4.1  If  f  g  Lip(l)  and  /Jn  tj  the  penodogram  from  a  realization  of  size  2  n  from  N(  0,  /), 
then 

1.  EYkn  =  log(ffj„  *  /(Jk/2n))  +  M/<?(logn/n)J 

2.  Eekn  =  MjO(log  n/n)3 

3.  Var Ykn  =  Var(f*n)  =  *3/6  +  A/30(iog  n/n)1 
I-  cov(Yhn,Y}n)  =  cov(ekn,e,n)  =  MjO( log  n/n)3 

uniformly  in  f  for  l<k<j<n  —  1. 

4.3  Properties  of  Wahba’s  estimate 

In  this  section  we  derive  properties  of  Wahba’s  estimate  which  allow  us  to  apply  the  results  of 
Chapter  3  when  we  use  this  estimate  of  the  spectrum,  in  other  words,  we  show  that  Wahba’s 
estimate  has  the  regularity  properties  which  Theorem  3.4  requires.  We  introduce  notation  and 
auxiliary  results  first. 

The  following  lemma  derives  properties  of  the  window  W*  m  „  and  its  derivatives. 

Lemma  4.2  If  Wxmn(a>)  is  given  by  (I-I),  0  <  k  <  2m,  and 

then 

s«P  KlnWl  = 

W  1  ' 

uniformly  in  A  >  0. 

Remark:  There  are  two  terms  involving  powers  of  A  in  the  expression  above  to  cover  A  — »  0  and 
A  — ♦  oo.  In  this  paper  we  will  usually  need  to  use  this  result  when  A  — *  0  in  which  case  Lemma 
4.2  implies 

Now,  notice  that,  for  «  large,  the  estimate  (4.3)  of  h  —  log  /  can  be  rewritten  as 

n 

hn(w)  =  ]T  (log(A-3„*/(fc/2n)+f*n)WA,m,n(w-k/2n) 

k=-n  + 1 
n 

=  £  HW(«  -  fc/2n)  log( Ajn  .  f(k/2n)) 

fc=-«+l 
n 

+  53  ~  */2 n)«*n 

fc  =  — n-f  1 

=  h,(w)  +  h;(w) 

so  that  it  consists  of  a  deterministic  term  and  a  random  term. 


(4.5) 
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4.3.1  Properties  of  the  deterministic  term 

In  this  section  we  find  bounds  for  the  deterministic  term  and  its  first  derivative.  In  addition, 
we  show  that  it  converges  in  Xj  norm  to  the  true  log-spectrum. 

Lemma  4.3  supw  |h„(w)|  <  16*1  +  (1  4-  A_1/,m)0(logn/n)  uniformly  in  A. 

Lemma  4.4  If  f  €  Lip(l)  is  differentiable,  then 

sup|h;HI  =  (l  +  A-l/Jm)0(l) 


uniformly  in  A. 

The  following  lemma  establishes  convergence  to  A  in  L j. 

Lemma  4.5  If  f  &  Lip(l)  and  h*  is  given  as  in  (4-5),  then 

||fc  -  Mi  <  4(A)  +  (1  +  A-l/,m)0(logn/n) 


uniformly  in  A,  i'-h.re 


A'(\) 


y*  /  M( 2ry)»"  \ 
V  1  +  A(2»l/),m  / 


1 


so  that  A(A)  is  bounded,  continuous,  and  A( 0)  =  0. 


(4.6) 

(4.7) 


4.3.2  Properties  of  the  random  term 

In  this  section  we  establish  bounds  for  the  random  term  and  its  first  derivative,  and  its  integrated 
mean  square.  First,  the  result  for  the  integrated  mean  square. 

Lemma  4.8  If  f  €  Lip(l)  and  h'  is  given  as  (4-5),  then 

E  j  h’Jwfdw  =  0(n-l)(l  +  A-l/,m). 

uniformly  in  A. 

The  following  result  establishes  the  regularity  properties  for  the  random  term. 

Lemma  4.7  Let  h*  (w)  =  5Z*=-«+i  f*»^n(<*')  and  let  the  following  conditions  hold: 

I.  Eckn  =  0(logJ  n/n3) 

S.  sap*£tof|*  =  0(1) 

S.  sup4ii  Efkntj*  =  0( log1  n/n3) 

4-  £kWn(u,-k/2n)=l 
5.  sup*  |W»(w)|  =  0(nr*~l) 

e.  mPjw:mi  =  o(«'*-1) 

7.  supw  |W;V)I  =  0(nr>-1), 

where  0  <  r,  <  1,  i  =  0, 1,  2.  Then,  for  any  0  <  6,s  <  1,  we  have  that 

s«P  KHI  <  Op(n1/3+<+'»"  +  n'+’’I_1)  (4.8) 

lot 
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This  result  requires  interpretation:  Assume  that  Wn  —  where  the  smoothing  sequence 

An  satisfies  A„  =  n~T .  Then  r0  =  r/2m,  rx  —  r/m ,  and  r 3  =  3r/2m.  Now  let 


/I  .  r  r  \ 

at  =  max  -  +  4  +  - - s,  s  + - 1 

\2  2m  m  ) 

and 

(\  .  r  3 r  N 

a3  -  max  -  +  <5  + - s,  s  + - 1  ■ 

\2  m  2m  J 

Then,  if  we  take  a  —  |  +  j  -  ^  and  0  <  r  <  m/3,  we  have  that  aj  =  1/4  +  3r/4m  +  <5/2  <  <5/2 
and  a3  =  — 1/4  +  5r/4m  +  4/2  <  i/2+1/6  so  that,  for  6  sufficiently  small,  it  follows  that  ax  <  0, 
a3  <  1/6,  and  hence,  sup|h^j  =  Or(  1),  and  sup|h*'|  =  0,(nl/,e). 


4.3.3  General  properties 

We  can  now  establish  a  result  for  the  integrated  mean  square  error  of  hn. 

Lemma  4.8  If  f  €  Lip(l)  and  h*  ia  given  by  (4-5),  then 

Ed,(h ,  kn)  =  A(A)  +  (1  +  A-,/,m)0(logn/n)  +  0(«'l/J)(l  +  A',''4m).  (4.10) 

uniformly  in  A. 

If  /«  =  exp(A„)  is  used  to  estimate  the  spectrum  /,  we  have 
Lemma  4.0  Let  f  €  Lip(l)  and  K  =  Kn~r  for  0  <  r  <  m/3.  Then 

dj(/n,  /)  =  0,(1)  (o(logn/n1"^)  +  A(tfn"')  +  0(n-l's)(l  +  A~,/4m))  . 
uniformly  in  A. 


4.4  Summary 

We  have  established  that,  if  /  €  Lip(l)  is  differentiable,  1 4* I  <  00 <  and  An  =  0(n_T), 
0  <  r  <  m/3,  the  estimate  of  the  spectral  density  given  by  /„  =  exp(h«)  has  the  following 
properties: 

1.  supw  |An(u>)|  =  0,(1)  (Lemmas  4.3  and  4.7)  which  implies 

(a)  there  exists  «  >  0  such  that  P( infw  /«(<*>)  >  e)  — *  1  and 

(b)  supw/„(w)  =  0,(1) 

2.  supw  |h'(<*»)|  =  0,(nl/e)  (Lemmas  4.4  and  4.7)  which  implies  supw  |/'(<*>)|  =  0,(n1/8) 

3.  ||/«  -  /||j  =  o(l)  (Lemma  4.9). 


These  are  sufficient  for  the  conditions  of  Theorem  3.4  to  obtain. 


Chapter  5 

Simulation  algorithms 


5.1  Introduction 

Recall  the  motivating  problem  stated  in  the  introduction:  a  realisation  X£~l  ~  jV(0,/),  / 
unknown,  is  observed,  procedure  0M  is  used  to  estimate  a  functional  0/  of  /,  and  it  is  of  interest 
to  estimate  the  sampling  distribution  of  0Z . 

In  chapter  3  we  proposed  the  following  approach  to  this  problem:  use  some  procedure  /n 
to  estimate  /  and,  by  simulating  many  realisations  from  a  process  {Yt}  whose  conditional 
distribution  given  /„  is  JV(0,/n),  determine  the  (conditional)  sampling  distribution  of  8V  and 
use  it  to  estimate  that  of  0X. 

From  a  practical  point  of  view  the  problem  is  two-fold:  (i)  finding  an  estimate  /„  of  /,  and 
(ii)  calculating  the  sampling  distribution  of  0y,  Y  ~  N(0,  /„).  The  first  problem  was  addressed 
in  the  previous  chapter  where  we  introduced  Wahba’s  (1980)  estimator  of  the  log-spectrum  and 
discussed  its  properties.  As  far  as  the  second  problem  is  concerned,  our  approach  will  be  to 
use  Monte-Carlo  to  approximate  the  distribution  of  0y  by  the  empirical  cdf  of  0y  for  many 
realisations  YJ*-1  ~  JV(0,  /). 

The  Monte-Carlo  approach  to  estimating  distributions  makes  generating  simulations  quickly 
very  desirable.  For  this  reason,  one  of  the  emphases  of  this  research  was  to  implement  fast 
algorithms  to  simulate  realisations  Y0n_l  ~  1V(0, g)  for  a  given  g. 

5.2  Simulating  Gaussian  stationary  processes 

There  are  several  ways  to  simulate  realisations  from  N( 0,  g).  For  example,  if  we  let  Co, . . . ,  cn-1 
denote  the  first  n  terms  in  the  Fourier  series  expansion  of  g,  then  (Lemma  2.1)  En  =  T(c0, . . . ,  c*_i) 
is  positive-definite  whenever  g  >  0.  It  is  thus  simple,  at  least  in  principle,  to  obtain  a  realisation 
from  Af(0,  E„):  let  Iq*-1  =  where  Z%~1  is  a  vector  of  iid  1V(0, 1)  variates. 

The  problem  with  this  approach  lies  in  calculating  Ej/3  Obtaining  the  square  root  of  an 
arbitrary  positive  definite  (p.d.)  matrix  requires  0(na)  operations.  Depending  on  the  amount 
of  structure  in  the  matrix  this  calculation  can  be  speeded  up.  For  example,  by  exploiting 
the  Toeplits  nature  of  En,  Carlin  et  al.  (1983)  describe  an  algorithm  (essentially  the  Levinson 
recursion)  to  obtain  a  Cholesky  factorisation  of  En  in  0(n3)  operations. 

The  question  is  whether  we  can  devise  an  algorithm  which  requires  less  than  0(n2)  operations 
to  generate  realisations  from  ff( 0,  g)  by  operating  in  the  frequency  domain  and  taking  advantage 
of  the  FFT’s  speed.  It  so  happens  that  there  is  an  algorithm  which  only  requires  0(n  log  n) 
operations.  The  constant  in  O  will  depend  on  the  smoothness  of  g. 
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5.2.1  Simulation  of  circular  processes 

In  order  to  explain  our  approach  it  is  easier  to  start  with  the  simplest  case.  It  turns  out  that 
this  is  the  case  where  the  covariance  matrix  is  a  circulant  £' ,  i.e,  {V,}  is  a  circular  process  with 
period  n.  We  saw  in  §2.4  that  E„  is  a  circulant  if,  and  only  if,  £'  =  F„A^F„,  where  F„  is  the 
Fourier  transform  matrix  and  A'  is  a  diagonal  matrix  whose  entries  are  the  DFT  of  the  first 
row  of  E'.  In  other  words,  the  singular  value  decomposition  of  a  circulant  can  be  implemented 
in  0(n  log  «)  time  by  taking  the  FFT  of  its  first  row. 

If  A'  is  p.d.  (i.e.  is  p.d.)  we  can  obtain  a  realization  yon_1  from  N( 0,£')  by  letting 
Y0n-1  =  Fny/\^Z^~l  for  Z*/1_l  ~  iV'(0,f),  Z0,  Zn/7  Ud  N(0, 1)  independently  of  the  other 
Z' s,  and  Zn^k  =  Zk,  where  for  convenience  we  assume  n  is  even.  Recall  that  Z\^~X  ~  JVc(0,  /) 
means  that  the  real  and  imaginary  parts  of  the  vector  1  ue  iid  N (0, 1/2). 

This  algorithm  can  be  implemented  in  0(n  log n)  operations:  (i)  O(nlogn)  operations  to 
obtain  the  FFT  of  the  first  row  of  E,  (ii)  O(n)  operations  to  calculate  its  square  root  and 
multiply  by  Zg~l,  and  (iii),  O(nlogn)  operations  to  obtain  the  inverse  FFT. 

In  implementing  this  algorithm  we  assumed  that  E'  is  a  covariance  matrix,  i.e.,  that  it  is 
not  negative  definite.  This  is  equivalent  to  assuming  that  the  DFT  of  the  first  row  of  E'  is  not 
negative  and  thus  each  of  its  elements  has  a  real  square  root. 

5.2.2  Simulation  of  non-circular  processes 

The  algorithm  used  in  §5.2.1  to  generate  realizations  from  circular  processes  in  0(logn/n) 
operations  can  also  be  used  to  generate  realizations  from  non-circular  processes  with  covariance 
matrix  E„. 

The  approach  in  this  case  is  to  generate  realizations  of  size  2n  from  a  circular  process  whose 
covariance  matrix  is  chosen  so  that  any  n-long  segment  of  the  2n-long  circular  process  has 
covariance  matrix  E„.  We  will  see  below  how  this  is  done.  This  idea  is  essentially  due  to 
Jonas  (1984)  except  that  in  our  implementation  we  introduce  a  correction  for  negative  definite 
matrices. 

Consider  the  2n  by  2n  Toeplitz  circulant 

=  T(co, ...  ,Cn _j,  0,  Cn_i,  ...  ,Ci).  (5.1) 

and  notice  that  it  has  the  property  that  any  n  by  n  block  originating  on  the  main  diagonal  is 
equal  to  En. 

Recall  that  since  E$n  is  a  circulant  it  has  the  factorization 

S',™  =  F3«AS^.  (5.2) 

In  other  words,  if  Aj„  is  p.d.  any  n-long  segment  of 

y0,n-1  =^A'Jnt/,Z0,"-^  (5.3) 

is  distributed  as  W(0, En),  where  Z?~l  ~  2Ve(0, 1),  Zo,  Zn  are  iid  N( 0, 1)  and  Z%  =  Zt for 
1  <  Jk  <  ». 

This  idea  will  only  work  if  E**  **  indeed  a  covariance  matrix,  i.e.,  if  it  is  not  negative  definite 
which  is  equivalent  to  the  entries  of  Aj„  being  non-negative.  In  fact,  whether  Ej„  is  p.d.  or  not 
is  a  function  of  the  smoothness  of  g  If  g  is  not  smooth  enough  it  may  happen  that  some  the 
entries  in  A|n  may  not  be  positive.  The  reason  for  this  is  easy  to  see:  Example  2  in  §2.4  showed 
the  entries  AJ  of  Ae  to  be  the  partial  sums  of  the  Fourier  series  of  g  at  the  Fourier  frequencies 
jfc/n,  k  =  0, . . .  ,2n-  1,  i.e.,  A*  =  *  j(k/n).  The  Dirichlet  kernel  is  not  positive  everywhere 

and  so  peaks  in  the  spectrum  may  be  weighted  negatively  and  leak  into  regions  of  low  spectrum 
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making  some  A£’s  negative.  On  the  other  hand,  if  g  €  Lip(l),  Lemma  2.20  ensures  that  there 
exists  some  N  such  that  for  every  m  >  N ,  DJm+l  *  g(X)  >  0  for  all  A. 

In  our  applications  g  =  /„  =  exp (hn)  >  0  where  An  is  Wahba’s  estimate  of  the  log  spectrum. 
Thus,  given  g  >  0,  we  can  find  Ejm  p.d.  by  determining  the  smallest  m  >  n  such  that  the 
eigenvalues  *  g(k/m)  >  0  for  fc  =  0 . m. 

In  practice  we  only  encountered  negative-definite  Ejn  when  simulating  short  realisations  of 
proceses  with  highly  peaked  spectra,  e.g.,  autoregressions  with  roots  near  the  unit  circle.  In 
order  to  avoid  too  many  increments  of  sise  1  in  searching  for  p.d.  E|m  we  adopted  the  strategy 
of  doubling  m  until  the  problem  disappeared.  In  other  words,  initially  we  simulated  series  of 
length  m  =  2*  where  the  actual  sample  sise  of  interest  is  n  <  2*,  i.e.,  we  constructed  for 
m  a  power  of  2.  We  then  calculate  the  DFT  of  the  first  row  of  £jm,  i.e.,  we  calculated  the 
eigenvalues  of  E|m.  If  some  turned  out  negative,  we  would  take  m  =  2*+l,  and  so  on,  until  all 
eigenvalues  were  positive. 


5.2.3  The  algorithm 

Recall  that  for  a  given  g  we  wish  to  assess  the  sampling  distribution  of  an  estimator  9y  of  9g  by 
Monte  Carlo.  We  detail  the  steps  of  the  algorithm  below. 

Notice  that  since  we  are  given  a  spectrum  g  and  not  covariances,  step  4  estimates  the 
covariances  by  approximating  the  integral  in  equation  (2.2)  by  a  sum. 

Specifically,  the  algorithm  consists  of  the  following  steps: 

1.  Get  g. 

2.  Set  d  •—  1  and  m  =  2*  >  n. 

3.  Set  d  —  2d. 

4.  Apply  the  DFT  to  g(k/dm),k  =  0, ...,  dm  -  1  and  obtain  c0, .... 

5.  Calculate  the  eigenvalues  AJ,  k  =  0, ....  dm  -  1  of  via  an  application  of  the  DFT  to 

(Co,  Cj  ,  ...,  —  l ,  0,  —  1 »  *-l  •) 

If  any  of  the  AJ’s  are  not  positive  go  back  to  step  2. 

6.  Simulate  a  realisation  Z  =  Z0, ....  £**-1  from  Nc(0,  A|m)  via  an  application  of  the  FFT  to 

*  =  0 . dm  -  1.  The  Uk  are  iid  IVC(0,1),  k  =  1 . dm/2  -  1,  Uk 

k  =  1  dm/2,  ...,dm  —  1,  and  U0,  Ujm/ j  are  simply  iid  N(0, 1). 

7.  Apply  the  FFT  to  Zo,  Zim.\  and  call  the  result  To.  •••,  Yim-i-  Because  Uk  = 
the  T’s  are  real  and  clearly  have  distribution  N(0,  Eeim). 

8.  Pick  *  such  that  0  <  »  <  dm  -  n,  and  thus  T0B_l  =  T,, ...,  Tn_i+,  is  a  realisation  from 

JV(o,En)  =  J\r(o,0). 

9.  Calculate  9y. 

10.  Go  back  to  step  6. 

After  repeating  steps  (6)-(10)  a  large  number  of  times  (100-1000)  we  usually  have  a  reasonable 
estimate  of  the  sampling  distribution  of  9y  under  g. 


Chapter  6 

Applications 


In  this  chapter  we  discuss  an  application  to  estimating  the  sampling  distribution  of  the  MLE 
(assuming  only  stationarity)  of  the  first  correlation  coefficient.  Specifically,  we  use  the  resampling 
algorithm  described  in  chapter  5  to  estimate  the  bias  and  variance  of  this  estimator.  We  then 
compare  the  estimates  obtained  by  resampling  with  Monte-Carlo  approximations  of  the  true 
finite-sample  bias  and  variance,  and  with  asymptotic  theory  approximations. 

It  is  important  to  note  that  the  asymptotic  approximations  require  knowledge  of  the  model 
and  for  this  reason:  (i)  will  not  be  very  useful  in  practice  and  (ii)  are  prone  to  be  more  accurate 
than  estimates  which  could  be  actually  used  in  practice;  they  are  presented  for  comparison 
purposes.  On  the  other  hand,  the  estimates  obtained  by  resampling  do  not  require  specific 
model  assumptions  beyond  stationarity,  Gaussianity,  and  spectral  density  in  Lip(l). 

We  present  results  from  simulations  under  two  different  models,  the  AR(1)  and  MA(1)  for 
several  parameter  values  and  sample  sixes. 


6.1  Estimating  correlation  in  the  AR(1)  and  MA(1)  mod¬ 
els 


Suppose  that  we  observe  a  realisation  X£~l  ~  W(0,  /)  and  we  wish  to  assess  the  sampling 
distribution  of  the  estimator 

-  _  /exp(»2*w)/*(w)du>  E. 

"  z,xi  (  1 

of 


f  exp(»2*w)/(w)dto>  /  cos(u»)/(to»)du> 

/  f(w)du  f  f(w)cL> 


under  two  different  models:  AR(1)  and  MA(1). 


6.1.1  AR(1)  Model 

The  AR(1)  model  is 

Xf  +  tfrXt-i  —  <i 

where  {<»}  are  iid  W(0,<rJ)  variates.  For  this  model  p  =  ^  and  (cf.  e.g.,  Kendall  and  Stuart, 
1976)  the  variance  of  the  limiting  distribution  of  y/n(p  -  p)  is 

V  =  (\-p>) 
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while  the  bias  of  p  is  given  by 

Ss-i+32. 

n  —  1 

8.1.2  MA(1)  Model 

The  MA(1)  model  is 

Xt  = 

where,  again,  {c(}  are  iid  N( 0,  <rJ)  variates.  For  this  model  p  =  ^/(l  +  V1*)  and  (cf.  e.g.,  Kendall 
and  Stuart,  1976)  the  variance  of  the  limiting  distribution  of  \/n(p  —  p)  is 

^  =  (1  -  3p5  +  4p4) 

while  the  bias  of  p  is  given  by 

B  ~  — ^(1  +  p)(4p*  -  2p  -  1). 

T1  —  1 

6.2  Results  from  simulations 

Below  we  present  results  Grom  simulations.  The  object  is  to  estimate  bias  and  variance  for  the 
estimator  of  the  first  correlation  given  by  (6.1).  We  compare  three  quantities: 

1.  asymptotic  theory  approximations 

2.  Monte-Carlo  approximations 

3.  bootstrap  estimates 

We  simulated  6  different  processes:  four  AR(1)  processes  with  parameters  .9,  -.9,  .5,  -.5,  and 
two  MA(1)  processes  with  parameter  -1  and  1.  For  each  of  these  models  we  considered  three 
sample  sixes:  32,  128,  and  512. 

6.2.1  Asymptotic  approximations 

The  asymptotic  approximations  were  calculated  using  the  formulae  in  the  previous  section.  We 
report  these  under  the  heading  oo. 

8.2.2  Monte-Carlo  approximations 

Monte-Carlo  approximations  of  bias  and  variability  were  obtained  for  each  model  and  each 
sample  sise  by  simulating  N\  —  10,000  realisations  of  the  process  and  reporting  the  average 
bias  and  average  square  deviation  from  the  mean  as  the  Monte-Carlo  estimates  of  bias  and 
variability.  In  other  words,  if  pi  denotes  the  estimate  of  the  value  of  the  first  correlation  p  based 
on  the  i**  sample,  we  estimate  bias  by 

1  S' 

f>Mc(p)  =  Pi  -  P=  P-  P- 

N 1  frt 

Similarly,  to  estimate  variance  we  used 

1  N' 

Va inc(p)  =  Yx  -  P)7- 

We  report  these  estimates  with  estimates  of  their  standard  errors  in  parentheses  under  the 
heading  MC. 
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6.2.3  Resampling  estimates 

Resampling  estimates  were  obtained  in  the  following  manner:  Ni  =  200  realisation  of  the 
original  process  were  simulated  and  for  each  an  estimate  of  the  spectrum  was  obtained  using 
Wahba’s  estimate.  In  addition,  for  each  of  these  estimates  of  the  spectrum  the  associated  first 
correlation  was  calculated,  i.e..  if  f\  is  the  estimate  of  the  spectrum  for  the  itk  realisation,  then 

.  _  Jc o$(2xu>)f\(u>)dw 
P'  ~  //;(-)(*- 


is  the  correlation  associated  with  /, . 

Then,  for  each  realization  (i.e.,  for  each  /’),  Ni  =  200  realizations  from  N(0,fn)  were 
simulated  and  we  estimated  bias  for  each  by 

bUc{p D  =  P'l  ~  p'  =  %  ~  P' 

and  variance  by 

Varvc(pD  =  £(#,  -  W)3 

where  p"}  is  the  estimate  of  correlation  from  the  jth  bootstrap  sample  of  the  i‘k  realization. 

Under  the  heading  BS  we  report  the  averages  of  these  quantities  over  the  200  realizations 
with  estimates  of  their  standard  errors  in  parentheses. 

6.2.4  Tables 

In  the  following  tables  we  use  the  notation  1.234_j  to  mean  1.234  x  10-J. 
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n  =  32 

BIAS 

VAR 

oo 

MC 

Is 

1.613., 

4.443.,  +  (0.152.,) 

3.392.,  + (0.162.,) 

2.344_, 

2.316.,  +  (0.042.,) 

2.612.,  +  (0.069.,) 

n  =  128 

oo 

MC 

BS 

0.394., 

1.212.,  +  (0.076.,) 
1.070.,  +  (0.048.,) 

0.586., 

0.584.,  +  (0.010.,) 
0.658.,  +  (0.015.,) 

EX3EQ3 

0.098., 

0.302.,  +  (0.038.,) 
0.307.,  +  (0.021.,) 

0.146., 

0.142.,  +  (0.0025.,) 
0.156.,  +(0.0024.,) 

Table  6.1:  Xt  +  O-SA^.i  =  et  =>  p,  —  —.5 


n  =  32 

BIAS 

VAR 

OO 

Me 

15 

1.613., 

3.054.,  +  (0.124.,) 

3.291.,  +  (0.116.,) 

1.563., 

1.525.,  +  (0.027.,) 

1.763.,  +  (0.048.,) 

n  =  128 

oo 

MC 

BS 

0.394., 

0.778.,  +  (0.062.,) 
0.816.,  +  (0.039.,) 

0.391., 

0.385.,  ±  (0.007.,) 
0.432.,  +  (0.009.,) 

n  =  512 

oo 

MC 

15 

0.098., 

0.178.,  +  (0.031.,) 
0.213.,  +  (0.017.,) 

0.098., 

0.094., +  (. 0017.,) 
0.109.,  +  (.0016.,) 

Table  6.2:  Xt  =  (t  —  1.0e»,  =>  pt  —  -.5 


n  =  32 

BIAS 

VAR 

oo 

MC 

BS 

-8.065., 

-4.291.,  +  (0.152.,) 
-3.520.,  +  (0.157.,) 

2.344., 

2.319.,  +  (0.043.,) 
2.552_,±  (0.072.,) 

n  =  128 

oo 

MC 

Is 

-1.969., 

-1.077.,  +  (0.077.,) 
-1.211.,  +  (0.053.,) 

0.586., 

0.591.,  +  (0.010.,) 
0.661.,  +  (0.016.,) 

BIB 

bum 

Table  6.3:  X , 
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Chapter  7 

Improving  estimators  by 
resampling 

7.1  Sufficiency  and  improved  estimators 

Suppose  that  we  observe  X£~l  with  distribution  P9  and  suppose  there  exists  a  sufficient  statistic 
Sn  =  5(Jlq_1)  for  ©.  In  this  case,  whenever  0  is  used  to  estimate  0,  9"  =  E&(9\Sn)  may  also 
be  used  to  estimate  9.  Notice  that  sufficiency  is  only  required  to  insure  that  9 *  does  not  depend 
on  0  and  is  hence  an  estimator. 

Now,  the  Rao-Blackwell  theorem  states  that  if  the  loss  function  L(9 , 0)  =  d(9  -  9)  is  convex, 
then  the  risk  function 

R(9",Q)  =  EqL(9‘,9)  <  R(9,&) 

for  all  0.  In  other  words,  if  there  exists  a  sufficient  statistic  S„,  C'  is  at  least  as  good  an 
estimator  as  9  and  possibly  better. 

For  example,  if  the  loss  is  quadratic,  i.e.,  L(9  -  0)  =  (0  —  0)3  the  risk  is  the  mean  square 
error: 

E9(9-9)2  =  Var(0)  +  b(0)3 

and  hence 

*(0*,e)  =  E9(E9(9\Sn)  -  9)*  =  R(9, 0)  -  £e(Var(0|S„)). 
so  that  the  reduction  in  risk  achieved  by  using  0*  instead  of  0  is 

£e(Var(0|Sn))  >  0.  (7.1) 

A  simple  argument  (cf.  e.g.  Lehmann,  1983)  shows  that  for  any  estimator  9(X*~l)  there 
exists  a  (possible  randomised)  estimator  9r  based  on  5„  with  the  same  risk  function.  Indeed, 
if  S„  is  sufficient,  the  conditional  distribution  of  X$~l  given  S„  does  not  depend  on  ©.  This 
implies  that  if  Y£~ 1  is  drawn  from  the  conditional  distribution  of  X£ -l  given  Sn,  the  risk 
function  of  (T  =  is  the  same  as  that  of  ^(XJ-1). 

The  estimator  S'  described  above  will  be  a  randomised  estimator  unless  it  is  equal  to  0  with 
probability  1.  In  this  case,  a  corollary  to  the  Rao-Blackwell  Theorem  asserts  that  0*  is  uniformly 
better  than  0 '  (and  hence  than  0). 

The  last  paragraph  motivates  the  approach  we  will  follow.  The  h*uristic  idea  is  to  find  a 
sufficient  statistic  Sn  of  lower  dimension  than  the  data  with  the  hope  that  9'  (and  hence  0) 
has  a  non-degenerate  distribution  given  S„  (and  hence  a  variance)  so  that  0*  is  indeed  a  better 
estimator  than  0  (compare  (7.1)). 
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7.2  Sufficiency  in  time  series:  the  circular  case 

If  the  data  are  a  realization  from  a  stationary  Gaussian  process,  i.e,  X ~  A7(0,En)  for 
Toeplitz  £n,  there  is,  in  general,  no  sufficient  statistic  of  dimension  lower  than  n  unless  we 
assume  the  data  are  generated  by  a  specific  model  with  a  reduced  number  of  parameters. 

The  types  of  models  commonly  used  for  time  series  data  assume  parametric  representations 
for  either,  (i)  the  mechanism  generating  the  data,  e.g.,  ARMA(p,q),  or  (ii),  for  the  spectrum.  The 
spectrum  in  general  does  not  determine  the  model,  but  it  does  under  the  additional  assumptions 
of  Gaussianity  and  stationarity. 

Our  strategy  for  finding  an  improved  estimator  8"  on  the  basis  of  8  is,  in  principle,  very 
simple:  apply  the  Rao- Blackwell  Theorem,  i.e.,  find  a  sufficient  statistic  Sn  under  the  assumed 
model  for  the  data  and  calculate  the  expectation  of  9  conditional  on  Sn.  The  problem  is  that 
under  most  of  these  models,  even  in  the  Gaussian  case,  it  is  not  easy  to  determine  a  reasonable 
sufficient  statistic,  let  alone  to  find  the  conditional  distribution  of  the  data  given  the  sufficient 
statistic. 

Because  of  the  difficulties  just  mentioned  we  will  take  a  different  approach:  we  will  exploit  the 
fact  that  the  problem  simplifies  enormously  under  the  assumption  that  the  process  is  circular. 
Under  this  assumption:  (i)  the  periodogram  is  a  sufficient  statistic  and  (ii),  the  distribution  of 
the  data  conditional  on  the  periodogtam  is  easy  to  obtain. 


7.2.1  Approximating  estimators  by  Monte-Carlo 

Assume  we  observe  XJ~l.  The  DFT  of  the  data  is  J*  =  F„X*~l  =  Ai/3r£~l 
diagonal  matrix  of  amplitudes  of  7*  and 

(  exp(»2iT7o)  \ 

rrl  = 

\  exp(»2*7n-i)  ) 


where  An  is  the 


(7.2) 


is  the  vector  of  phases.  Since  XJ-1  is  real  we  have  that  A;  =  An_,  and  -r}  =  -yn_,. 

Under  the  assumption  of  circularity,  y,  are  iid  17(0, 1)  for  j  =  1, . . . ,  n/2  —  1.  This  means 
that,  given  the  periodogram,  the  amplitude  of  the  DFT  is  fixed  and  all  we  have  to  do  to  obtain 
a  sample  yon_1  with  the  same  conditional  (on  the  periodogram)  distribution  as  the  original 
data  Xq  ~l  is  to  draw  iid  17(0, 1)  phases  Aj-1,  multiply  each  exp(t2ir£,)  by  the  corresponding 
amplitude  A^ J  (given  by  the  square  root  of  the  periodogram),  and  finally  take  the  inverse  DFT. 
Using  this  device,  many  new  samples  with  the  same  conditional  distribution  as  the  original  data 
may  be  generated  and,  for  each,  the  value  of  8  can  be  calculated  and  averaged  with  the  rest  to 
yield  an  arbitrarily  precise  approximation  9‘  to  EQ(9\Sn). 

The  obvious  question  which  arises  is:  why  go  to  all  this  trouble?  Why  not  simply  start 
off  with  an  efficient  estimator  (e.g.,  the  MLE)  and  avoid  (usually  expensive)  resampling?  The 
answer  is  that  efficient  estimators  may  not  always  be  known,  or  may  not  have  simple  represen¬ 
tations,  but  there  may  be  other  estimators,  albeit  inefficient,  which  are  easy  to  calculate.  This 
point  is  illustrated  with  an  example. 


7.2.2  Example:  crossing  rates 

Suppose  we  are  interested  in  estimating  the  rate  <j>a  at  which  a  process  crosses  a  barrier  at  height 
a,  i.e.,  the  rate  at  which  the  process  =  X*  -  a  changes  sign.  In  other  words,  if  we  define 


1  if  Xk  >  a 
0  otherwise, 


(7.3) 
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we  wish  to  assess 

1  n"’ 

<t>a  =  lim  - -  i  -  Zi,)2- 

n  —  oo  n  —  1  ' 

k  =0 

Since  (Zk  + 1  —  A*)3  =  1  if,  and  only  if,  a  crossing  of  a  ocurrs  at  time  k  +  1,  the  problem  reduces 
to  evaluating 

<j>a  =  E(Zk+ 1  -  Zk)2 

=  EZ2i+l+EZl-2EZk+lZk 
=  2P(A*  >  a)  -  2P(Xk+i  >onX»  >a) 

=  2*(a/<rz)  -  2P(Xk+l  >  a  n  Xk  >  a) 

Even  in  the  Gaussian  case  there  is  no  simple  closed  form  representation  for  <pa.  However,  in 
the  special  case  when  a  =  0,  we  have  that  (cf.  e.g.  Johnson  fc  Rots,  1972) 

2P(Xk+i  >  0  n  X*  >  0)  =  .5  +  *~l  arcsin(p) 

where  p  —  EX0Xi/ EXg,  so  that 

4>o  =  .5  —  r_l  arcsin(p). 

An  efficient  estimator  of  p  is  the  MLE  (under  the  Toeplits  and  the  AR(1)  models)  p  = 
Y  XkXk+i/  Y  Xk  so  that  <t>o  =  .5  -  t-1  arcsin(p)  is  the  MLE  of  <j> o- 
On  the  other  hand,  consider  the  naive  estimator  of  <j>a  given  by 

0 

Notice  that  <j>"a  can  be  calculated  easily  for  all  a  but  it  is  clearly  inefficient  since  it  uses  only  the 
information  in  the  sign  of  {A*  -  a}. 

As  we  mentioned  above,  under  the  assumption  that  our  process  is  circular  we  can  resample 
the  phases  of  the  DFT  of  our  original  realisation  to  construct  an  arbitrarily  close  approximation 
to  a  third  estimator  of  <fra, 

=  e9  (r\r*) 

where  /*  is  the  periodogram  of  the  data  X$~l  evaluated  at  the  Fourier  frequencies. 

Two  questions  arise  now:  first,  what  are  the  relative  efficiencies  of  these  three  estimators 
and,  second,  in  what  situations,  if  any,  can  we  improve  our  estimators  by  resampling  if  the 
original  process  is  not  circular?  We  will  defer  addressing  the  first  question  until  the  simulations 
and  try  to  shed  some  light  on  the  second  one. 

7.3  Improved  estimators:  non-circular  case 

Our  approach  is  suggested  by  that  of  the  circular  case:  calculate  the  periodogram  of  the  original 
data  and  estimate  the  random  variable  E(6\I*)  by  8’  where  8"  is  obtained  by  resampling  our 
original  data  conditional  on  the  periodogram  and  pretending  that  the  phases  are  iid  17(0, 1)  to 
obtain  new  samples.  For  each  new  sample  we  then  calculate  the  value  of  8  and,  finally,  take  the 
mean  of  the  values  of  8  for  each  sample  to  be  our  approximation  to  8‘ . 

Since  in  the  non-circular  case  the  phases  are  not  really  iid  17(0,1),  but  only  approximately 
so,  we  will  have  that  8 *  ^  E{8\I*}  in  general  (the  r.h.s.  may  not  even  be  a  statistic  since  it 
may  depend  on  8),  but  if  the  difference  is  not  large  we  may  still  get  an  improvement  over  8.  For 
example,  if  8  is  unbiased,  then  8‘  will  in  general  be  biased  but  will  have  lower  variance  so  that 
in  the  end  the  MSE  may  be  reduced  by  using  8 *. 
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7.3.1  Distribution  of  resampled  data 

In  order  to  study  more  precisely  what  effect  pretending  the  phases  are  iid  1/(0, 1)  has  on  the 
resampled  data  V0n_l,  we  determine  the  distribution  of  yon  ~ 1 . 

Lemma  7.1  Let  X$~l  ~  lV<0,En)  ,  En  =  T(co . e«_i),  and  J*  =  F„Xq~1  = 

where  An  is  the  diagonal  matrix  of  amplitudes  I*(k/2n),  k  =  0, . . .  ,n  —  1.  If 

(exp(t2*A0)  \ 

exp(t2s,6„_i )  / 

where  the  components  of  1  are  ltd  (7(0,1),  Si,  —  bn-fc  for  0  <  k  <  n,  and  FJ*  1  = 

T^aJ/’a;;-1,  then 

C(Y0'-l\Z)  =  C(XZ~l\X  =  Z) 

where  X£-1  ~  N( 0,£‘),  =  T(c$ . and  c‘  =  (1  -  !*|/n)c*  +  (|fc|/n)cn_|t|  for 

0  <  *  <  n  -  1. 

Remark:  Lemma  7.1  establishes  that,  for  n  large,  the  covariances  c£  of  the  data  resampled, 
pretending  the  process  to  be  circular,  approximate  those  of  the  original  data  for  k  not  too  large. 


7.3.2  Heuristics 

We  can  now  study  what  9’  estimates  in  the  non-circular  case. 

Assume  that  the  parameter  of  interest  is  a  function  of  the  unkown  covariance  structure 
(co,  ci, . . .)  of  the  process,  i.e.,  9zn  =  9(f).  The  corresponding  parameter  for  the  circular  process 
is  9/n  =  9(f„)  where  f„  denotes  the  following  “sum”  of  the  first  n  Fourier  coefficients  of  /,  i.e., 

n  -l 

/«(*)=  Yh  ((!  -  lfel/n)c»  +  (lfc!/n)c"-l*l)elP(,2lru’fc)- 

fc=—  n+1 

Assume,  for  simplicity,  that  9  is  unbiased,  i.e.,  E(9)  =  9j,  so  that  its  MSE  is  just  its  variance, 

E(9-9,)2  =  Var(0). 

On  the  other  hand  the  distribution  of  the  resampled  data  is  N(0,  /„)  and  since  9  is  unbiased 
it  follows  that  E(9 *)  —  9fm,  and  hence  the  MSE  of  9 *  is 

E(6'  -9, )'  =  V*t(9-)  +  (9/n  -9,)\ 

It  remains  to  assess  bn  =  9f„  —  9f.  To  do  this  it  will  be  convenient  to  consider  9  as  a  function 
of  the  Fourier  coefficients  of  /,  i.e., 

9(f)  =  0(co,ci, ...). 


Similarly 

0(/n)  =  0(co,(l  -  l/n)ci  4-  Cn_i/n, . . . ,  (1  -  l/n)ci  +  c„_!/n,  0, . . .), 
so  that,  using  a  Taylor  expansion,  it  is  easy  to  see  that 


bn 


~  6 1 
n-l 


-E 


*=o 


!*!>» 


09 
Ch  dck 
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The  quantity  bn  measures  the  change  that  circularising  by  assuming  iid  random  phases 
induces  in  the  parameter  of  interest.  If  )6n|  is  not  too  large  and  8  is  rather  inefficient,  we  may 
do  substantially  better  by  using  9’  rather  than  9  even  in  the  non-circular  case. 

In  the  example  concerning  the  estimation  of  the  crossing  rate  of  0,  <p0  =  -5  -  »_l  sin_1(p), 
p  =  co/ci,  and  it  is  straightforward  to  verify  that 


n*s/l  -  p2 

and  the  results  from  our  simulations  indicate  that  there  is  considerable  reduction  of  MSE.  In 
fact,  it  turns  out  that  is  almost  as  efficient  as  the  MLE  even  in  cases  where  the  parameter 
space  is  of  much  smaller  dimension  than  the  periodogram.  The  reason  for  this  is  that  the  bias 
introduced  by  assuming  circularity  is  of  the  same  order  of  magnitude  as  the  bias  of  the  MLE 
and  thus  contributes  little  to  the  MSE  for  n  large. 


7.4  Results  from  simulations 

In  this  section  we  report  results  from  the  simulations.  For  each  process  we  simulated  N\  =  1,  000 
realisations  and  estimated  phio  by  fa,  4> Oi  and  fa,  reported  under  the  headings  naive,  MLE,  and 
resampled,  respectively.  The  resampled  estimator,  fa,  was  obtained  by  resampling  the  phases 
of  each  realisations  100  times,  i.e.,  it  is  the  average  of  100  fa's. 


» 


CHAPTER  7.  IMPROVING  ESTIMATORS  BY  RESAMPLING 


40 


n  =  512 

BIAS 

MSE 

Naive 

MLE 

Resampled 

bn 

-0 

-0.360.,  +  (.048.,) 
-0.184.,  +  (. 047.,) 

-0.128_i 

0.520.3  ±(.030_3) 
0.225.3  ±  (.013_3) 
0.225.3  ±(.013_s) 

n  =  128 

-0 

-1.141.,  ±(  098_i) 
-0.464.,  +  (.099.,) 
-0.513., 

2.022-s  ±(.114_s) 
0.962.3  +  (  063.3) 
0.975_,  +  (.061_3) 

Naive 

MLE 

Resampled 

fcn 

Tabic  7.1:  Xi  4- . 9 ^  =>  <£q  =  .8564 


n  =  128 

BIAS 

MSE 

Naive 

MLE 

Resampled 

bn 

0 

1.066.,  +  (.095.,) 
0.349-,  +(. 097.,) 
0.513_, 

2.093.3  +  (.112.$) 
0.893.,  ±  (.055.*) 
0.945., +  (. 056.,) 

Tabic  7.2:  Xt  —  Ci  ^  <f>Q  ~  .1436 


Appendix  A 

Auxiliary  results  and  proofs 


A.l  Chapter  2 

Proof  of  Lemma  2.1:  Fix  n  and  let  =  T(co, . . . ,  c,i_i).  For  any  vector  (*o>  ■  •  • ,  xn- 1)  define 
<^X(A)  =  *•  exp(— «2tsA).  Then,  since  /  >  0  and  <j>  is  continuous, 


0  <  ['  |rf.(A)|a/(A)«tt  =  E  =  *'£ 

■'°  »,i=o 


unless  |<0*(A)|  =  0«l,  =  0. 
Proof  of  lemma  2.2 


l/^MI  =  I  E  c«(*2,r*)*  exp(»2x»u<)|  <  (2ir )*  E  I 

which  proves  the  first  implication  of  the  first  part  of  the  lemma.  As  for  the  second  implication, 
it  is  easy  to  see  that  if  a  function  /  has  a  bounded  first  derivative  then  Mf  =  sup,  |/'(*)|  and 
hence  /  6  Lip(l). 

To  show  /  continuous:  Let  t  >  0.  Since  £|c,|  <  00  there  exists  N  =  N(f)  such  that 
E|,j>AT  lc*l  <  f/2-  80  that 


|/(«  +  a)-  f(v)\  = 


< 

< 


E^.O  -  exp(i2x*5))  exp(»2xsu>) 

t 

2  E  lc*  +  « 

l»l<w 

2x6  E  MM  +  « 

■<» 


but  for  fixed  AT,  the  first  term  above  is  continuous  in  6  and  goes  to  0  as  6  — *  0.  ■ 

Proof  of  lemma  2.4:  The  first  assertion  follows  immediately  from  the  properties  of  Z(w).  As 
for  the  second,  notice 

\F  -  G\*  =  f  +  g  -  2y/fgcos(6  -  <f>) 

where  8  and  4>  are  the  phases  of  F  and  G  respectively.  Clearly  9  =  <p  achieves  the  minimum.  ■ 
Proof  of  Lemma  2.5:  Obvious.  ■ 
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Proof  of  Lemma  2.10:  The  first  assertion  is  standard  (cf.  e.g.  Hasminskii  and  Ibragimov, 
1986).  As  for  the  second,  let 

An  =  jy"{/(z)  -  f(*)}Dn(t  -  t)Dn(x  -  »)dl 

=  j  j {/(z  A  t)  -  /(s)}D«(z  +  t  -  i)Dn(z)dz 

<  j 1 f(x  At)-  f(s)\\Dn(x  A  t  -  s)||Dn(x)|dx 

<  Mf  J  min(|x  +  t  -  s|,  |1  -  *  -  t  +  s|)|Dn(z  +  t  -  s)||D„(x)|dz 

since  /  €  Lip(l),  periodic  and  symmetric,  so  that 

An  <  Mf  sup  |zDn(x)|  [  j Dn(x)\dx. 

0<z<t  J 

Since  f  \D„(x)\dz  =  O(logn)  (cf.  e.g.,  Butler  and  Nessel,  1971),  to  establish  the  lemma  it 
suffices  to  show  that  there  exists  K  such  that  sup0<z<  t  \xDn(x)\  <  K.  Indeed, 


sup  \xD„(x)\  =  sup 


|  z  sin(mx)  I 


o<*<.s  |  sin(irx)  | 


<  sup 


0<r<-5  |  sin(xz)  I  5x 


since  (cf.  Butser  and  Nessel,  1971), 


sin(xz)  1  -  x1 
xx  ~  1  +  x3 


and  the  lemma  follows  immediately.  ■ 

Proof  of  lemma  2.11:  The  first  assertion  is  obvious.  The  second  is  a  standard  result  (cf.  e.g. 
Priestley,  1981)  and  is  also  a  consequence  of  the  proof  of  Lemma  2.12  which  follows.  ■ 

Proof  of  lemma  2.12:  Let  A  =  l/n  and  M  =  m/n.  Using  the  Cramer  representation  for  Xt 
and  the  fact  that  £)"_T0l  exp(i2xAt)  =  exp{i*(n  -  l)A}Dn(A)  it  is  easy  to  see  that 

j;M  =^J  F(w)€*'<B-1><“-''>£>n(w  -  P)dZ(v).  ( A- 1 ) 

Using  the  fact  that  F(u>)  =  F(-u>),  the  symmetry  of  Dn,  and  straightforward  manipulation,  we 
have  that 

»{';(**)}  =  \{K{p)  +  w>] 

=  2^/n  /  F(“')IexP(',r(n  -  W"  ~  P))D*(W  ~  M) 

+  exp(ix(n  -  l)(w  +  p))Dn(u>  A  n)]dZ{v) 


=  ^{JliP)-JnM} 


~  2^n  f  f(u'^exp("r(n  ~  *)("  -M))0„(«  -m) 

—  exp(ix(n  -!)(<*  +  fi))Dn(v  A  n)]dZ(ui). 
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Using  (2.3)  and  the  symmetry  of  /  and  Dn,  it  is  straightforward  to  see  that 

a,  =  a)} 

=  ^  J  /(“'){« P(**(n  -  1)(— ^  +  A))Z>„(w  -  fi)Dn(u>  -  A) 

+  exp(iir(n  -  l)(-p  -  A))Z?n(w  -  p)D„(w  +  A) 

+  exp(i»(n  -  A))Un(u<  +/i)Dn(w  -  A)) 

+  exp(iir(n  -  1)(m  -  A ))£>„(<*>  +  A)}du> 

-  “M*  -  *»  /  /mp.(,  -  »)P.(»  -  »*, 

=  y  f(v)Dn(ui  -  -  A)t. 

y  /hm* + -  a)j», 

and  the  result  follows  from  lemma  2.10. 

Similarly, 

fln  =  E^J^tSJ^X) 

=  y  /H{exp(»ir(»  -  l)(-/4  +  A ))Dn(*  -  /t)£.(w  -  A) 

-exp(»ir(n  -  l)(-p  -  A))D„(u>  -  fi)Dn(w  +  A) 

-  exp(»w(n  -  1  )(p  +  A ) ) Dn (u;  +  n)D„{u>  -  A)) 

+  exp(»T(n  -  1)(**  -  A ))Dn{u/  +  t*)Dn(u  +  A)}du; 

-  «■»(;(■»  -  >)(»  -  •»)  y /(a)D.(.  _  ,)D.(„  _  a),^. 

1)(,A  A))  y/MD>(„+„)D>(„_A)4, 

=  A>)  y  -  A)4. 

_co.(«ua-a))  y  + 

and 

Cn  =  EXJZMSJX  A) 

=  J  /(u/){exp(iT(n  -  1)(— /a  +  A))D„(u>  -  -  A) 

-  exp(»*(n  -  l)(-p  -  A))Z?n(ui  -  +  A) 

+  exp(«T(»  -  l)(/a  +  A ))Dn(u>  +  n)Dn(u  -  A)) 

-  exp(i»(n  -  l)(#s  -  A))D„(u>  +  n)Dn{w  +  A)}du> 

- - i - — -  /  f(w)Dn(w  -  m)D„(w  -  \)du> 

sin(*(n  -  l)(u  +  A))  f  „  ,  _  .  .  _  .  ... 

■+ - — - J  f(w)Dn(v  4-  n)D„(w  —  A)dw 
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=  —*$n  -A))  /  /(«)*>-(«  -  -  A)*- 

~S1°~  2n+  A))  /  /(«)*>•(«  +  -  A)^ 

=  MjO(log  n/n). 


Proof  of  lemma  2.13:  To  prove  the  first  assertion,  let 

Vi  =  -^={A’v/r+7  +  y'v/irrp} 

and 

Ui  =  -L  |jf  jr~P  -  r  s/TTp) 

for  X,Y  iid  JV(0, 1).  It  is  easy  to  verify  that  Ui,Uj  are  also  iid  #(0, 1), 

X  =  -j=  [uls/\  -r  p  +  Ui\Jl  —  p| , 

Y  =  — 7=  {tfj  y/\  —  p  —  Ui\/\  +  p| , 

and  W  —  (1  +  p)U*  +  (1  -  p)Uj .  This  proves  the  first  assertion. 

To  prove  the  second  assertion  it  suffices  to  show  that  the  random  variables  Z\  —  U*  +  U\ 
and  Zi  —  (U*  -  U\)l{Ul  +  U\)  are  independent. 

Indeed,  verify  U\  —  Zx(  1  —  Zj)/2,  U*  =  Zi(  1  +  Zj)/2.  Then,  verify  the  Jacobian  is  Zi/2. 
At  last,  verify  the  resulting  joint  density  of  Z\  and  Z\  is  given  by 

ir-1exp(-*t/2)(l-*5)-l/* 

and  the  result  follows.  ■ 

Proof  of  lemma  2.14:  Using  the  Cramer  representation  (2.3),  it  is  straightforward  to  see  that 

*D,(«o  -  A)D»(w0  ~  X)Dn(<t>i  -  p)Dn(ij»l  -  p)dZ(u>o)dZ(w'0)dZ(u>i)dZ(u[). 

Taking  expectation  yields  contributions  in  three  cases: 

1.  u>o  —  w0  and  ~ 

2.  wo  —  Wj  and  wj  =  w0 

3.  u»o  —  — <*»i  and  Wq  =  — Wj. 

The  contribution  from  (1)  is 

±/ f(v0)Dl(w0  -  A)du>0^  J  g(u>i)Dl(wi  -  p)du/i  = 

f  f(u>)Kn(v  -  X)dw  J  g(w)Kn{v  -  X)du>  =  EIl(X)EPn(p). 

The  contribution  from  (2)  is 


\_ 

n 


J  F(w0)G(wo)e,’<B-1><<'-MDn(u;o  -  A)D„(w0  -  p)dw0x 
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xiy*  G(Wl)F(Wl)e-,>*-1«'“J1Dn(«1  -  A)D.(Wl  = 

1  I  f  -  J 

=  \J  F(w)G(w)Dn(u  -  h)Dn(ui  -  fi)dv  . 

Analogously,  the  contribution  from  (3)  is 

J  F(v)G(Z)Dn(u-\)Dn(u,  +  »)dv  . 

The  result  now  follows  trivially  ■ 

Proof  of  lemma  2.16:  Consider  each  of  the  terms  in  (2.9).  By  Lemma  2.10,  the  first  term 

J  f(x)Dn(x-X)Dn(x-n)dx  =  f(n)Dn(,i*  —  A) 

+  J  {/(*)  '  ~  >)£»(*  -  H)dX- 

and  a  similar  expression  holds  for  the  second  term.  Ifp  ^  A  are  Fourier  frequencies,  0  <  ft,  A  <  .5, 
D„(fi  —  A)  =  Dn(n  +  A)  =  0  and  Lemma  2.10  gives  the  desired  result.  ■ 

Before  we  can  prove  lemma  2.17  we  need  three  auxiliary  results.  First  a  result  bounding  the 
distribution  of  quadratic  forms: 

Lemma  A.l  If  X  ~  AT(/i,  E),E  >  0,  then,  there  exists  a  random  variable  W  ~  x\ 
such  that 

t.W  <  X'X  <  i"W 

where  t.,tm  are  the  smallest  and  largest  eigenvalues  ofE  respectively. 

Proof:  Let  E  =  JTJ'  where  T  is  the  diagonal  matrix  whose  entries  are  the  eigenvalues  of 
E  and  J  is  the  matrix  of  normalised  eigenvectors  such  that  J  J'  —  J'J  =  I.  Notice  that 
X  =  JTl/2(T~l'2J'n  +  Y)  for  some  Y  ~  1V(0,  /)  so  that 

X'X  =  (Y1  +  n'JT-lf2)T(T-l/tJ'n  +  Y)  =  Z'TZ 

for  Z  =  T-Wj'p  +  Y  ».  N(T~x,2J'p,  I). 

But  Z'TZ  =  ^iiZf  so  that  clearly 

t.  ]T  Z2  <  Z'TZ  z> 

and  the  result  follows  from  the  fact  that  ^  Z2  ~  xi(E(Z)' E(Z))  and  that  E(Z)'E(Z)  = 

Second,  a  lemma  which  bounds  the  expectation  of  the  log  of  a  non-central  x\- 

Lemma  A. 2  If  W  ~  Xj(203)  is  a  random  variable  with  ron-central  x2  distribution  with  2 
degrees  of  freedom  and  non-centrality  parameter  202 ,  then 

-7<  Elog(W72)<03(l+T)-7 

where  7  =  .57721 ...  is  Euler’s  constant. 

Proof:  Biekel  and  Doksum  (1977)  show  that  a  x\(202)  variate  can  be  generated  by  first  sampling 
a  Poisson  (62)  variate  R  and  then  sampling  a  central  *5«+j  variate  independently.  Also,  it  is 
well  known  that 


£l°g(Xj,/2)  =  =  ib(r) 
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where  is  the  digamma  function  (cf.  e.g.  Abramowits  and  Stegun,  1954).  In  particular, 
V>(1)  =  -7,  so  that 


E( log  W/ 2) 


oo 


E 


exp(-g3)03V(lfe-H) 

fc! 


exp(— <9a)(V>(l)  +  ^ 

k  =  \ 


»u4>(k  +  l) 

k\ 


Since  0  <  tp(k  +  1)  <  k,  for  k  >  1  (cf.  Abramowit*  and  Stegun,  1954),  all  the  terms  in  the 
sum  are  positive  so  that  the  sum  itself  is  positive  which  implies 

exp(-03)tf(l)  <  £log  W/2. 

Since  V(l)  =  -7  and  0  <  exp(  —  93)  <  1  the  first  inequality  follows. 

The  second  follows  from  the  same  property,  0  <  i>(k  +  1)  <  k,  k  >  1,  so  that 


OO 


E 


eu^k  +  i) 

k\ 


»  nih 

<  y - 

1)! 


03exp(03). 


Thus,  we  have  that 

£logW/2<  -7exp(-fl3)  +81  <  -7(1  -  63)  +  93  =  -7  +  02(l  +7)-« 

Lemma  A. 3  If  U  and  V  are  two-dimensional  vectors  such  that 

(  V  )  ~*(°-(  rA'  Tl  ))' 

where  r3  is  the  largest  eigenvalue  of  r3  A!  A  and  h(U)  is  a  real- valued  random  variable  with  finite 
second  moment  er3 ,  then  cov(h(U),  log  V'V)  =  trO(r3),  uniformly  in  A  and  h. 

Proof  of  Lemma  A. 3: 


cov(h(U),logV'V)  =  cov(h({7), log  V'V/2)  = 


E{(h(U)  -  Eh(U))( log  V'V/2  -  E log  V'V/2 )}  = 

E{{h(U)  -  Eh(U))E(\os  V'V/2  -  Elog  V'V/2)\U}  .  (A.2) 

First  we  evaluate  the  conditional  expectation  above.  We  know  E( log  V'V/2)  =  -7  and  the 
distribution  of  V  given  U  is  simply  N(tA'U,I  -  r3A'A).  Write  A' A  =  PAP',  for  orthogonal  P 
and  diagonal  A  with  entries  1,  A],  where  0  <  Aj  <  1. 

First  we  bound  V'V\U\  apply  lemma  A.l  with  E  =  /  —  r3  A'  A  =  P(I  -  t3A)P'  and  /s  = 
rA'U  =  rPAl^3QU ,  where  Q  is  orthogonal.  If  we  let  QU  =  R,  it  is  easy  to  verify  that  the 
non-centrality  parameter  is  given  by 


203  =  t3U'Q'  A1/3(I  —  t3 A)~l A^^QU 


<  R'R- 


=  U'U- 


1  -  t3  1  -  r3 

Since  the  eigenvalues  of  I  —  r3A'A  are  given  by  1  -  r3  =  t.  <  t‘  —  1  —  t3Aj  it  follows  that 

(1  -  r3)W  <  V'V\U  <  (1  -  r3A,)W,  (A. 3) 


APPENDIX  A.  AUXILIARY  RESULTS  AND  PROOFS 


47 


for  some  W  ~  *3(203).  Taking  logs  and  expectation  in  A. 3,  an  application  of  lemma  A. 2  yields 
log{l  -  T3)  <  E  u'j  +y  <  log(l  -  r3A2)-f-  ^  -  ^(1  +7) 

which  implies 

«  u)  +T|  <  >o8  (j-hj)  ♦  +  „  .  +T)  (H2  +  ,)  o„., 

Lemma  A. 3  follows  upon  using  this  result  and  applying  Schwars's  inequality  to  evaluate  (A. 2) 
together  with  the  fact  that  E(U'U)1  =  Ex j  =  8  and  h(U)  has  standard  deviation  cr.  ■ 

We  can  now  prove  lemma  2.17:  By  lemma  2.13  the  two  periodogram  ordinates  can  be  written 
as 


Wi  =  (tf1s  +  tr3)(i  +  PlS1) 
=  w-(  l  +  MO 


where  p i  =  MjO(logn/n),  and 


Since  |0,|  <  1 


Wi  =  (V>+V’)(l  +p202) 
=  ^'(1  +Pj«S). 


log(W.)  =  log(W,*)  +  log(l  +  px0x) 
=  log(W,-)+p,0,+O(p3). 


so  that, 


cov(log(Wj),log(Wj))  =  cov(log(Wi),log(W;))  +  p,cov(log  W’,82) 

+Picov(log  +  0{\pi\  +  |pi|)3. 


(A.4) 


By  lemma  2.12  the  last  term  above  is  Af30(log  n/n)3. 

It  remains  to  calculate  each  of  the  covariances  in  (A.4).  Notice  that  U\  and  Ut  axe  in¬ 
dependent  and  so  are  Vj  and  Vj,  so  that  the  covariance  matrix  of  U\,  C7j,  V),  and  V3  looks 
like 

U  7)- 

Lemma  2.12  tells  us  that  r,  plt  and  p3  are  MjO(\o%n/n)  and  an  application  of  Lemma  A. 3 
yields 

cov(log(Wt),  log(Wj))  =  0(r3)  =  MjO{  logn/n)3.B 
Proof  of  lemma  2.18:  Clearly 


FBEF„  =  (-  V  <rk,txp(i2*(tl  -  jk)/n) 

\nk,.=o  j 


(A.5) 


Equation  (2.13)  follows  easily  using  l  =  j  in  (A.5)  and  the  fact  that  <rs»  =  c\k-,\-  Equation 
(2.14)  follows  upon  dividing  (A.5)  into  two  sums,  £o~'  and  Sln+i'  changing  the  summation 
index  in  the  second  to  t'  =  n  +  t,  and  using  tbe  fact  that  <7*  =  <r_*.  ■ 

Proof  of  lemma  2.19:  Immediate  from  the  fact  that  E,n  is  a  circulant.  ■ 
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Proof  of  lemma  2.20:  By  Lemma  2.2  /  is  continous  and  since  it  is  positive  /(w)  >  e  for  some 
e  >  0.  Now, 


/(w)  =  c*  exp(i2*wifc)  +  ^  ck  exp(.2irwfc) 

lk|<n  |k|>» 

=  an(ui)  +  rn(u>) 

Notice  that  an(k /2n)  =  Xck  is  the  kth  eigenvalue  of  and  that  a„  and  rn  are  real- valued 
functions.  On  the  other  hand,  since  ^  |c*|  <  oo  there  exists  N  =  N(e)  such  that  rn(w)  <  e/2 
for  n  >  N  and  thus 

On(u/)  =  f(u> )  -  rn(w)  >  e  -  e/2  =  e/2. 

which  means  that  all  the  eigenvalues  of  Ejn  are  positive  for  n  >  N .  ■ 

A. 2  Chapter  3 

A. 2.1  Auxiliary  definitions  and  results 

This  section  is  devoted  to  establishing  definitions  and  auxiliary  lemmas  which  we  use  in  deriving 
the  results  in  Chapter  3. 

In  what  follows,  let  fi  be  complete  and  separable.  Let  B  be  the  Borel  cr-algebra  and  M  the 
space  of  probability  measures  on  (0,  B).  Let  d(x,y)  =  ||x  -  y||  denote  the  metric  on  fl. 

Definition  A.l  (Huber,  1981)  For  any  subset  A  C  0  define  the  closed  S -neighborhood  of  A  as 
A‘  =  {x  €  fi  I  inf,eA  d(x,  y)  <  6}. 

Definition  A. 2  (Huber,  1981)  The  Prohorov  distance  between  two  members  F,G  €  Ad  ts 

dP(F,G)  =  inf  {«  >  0|F{A}  <  G{A*} +  c  for  all  A  €  B}  . 

Theorem  A. 4  (Strassen)  (Huber,  1981)  Let  e  >  0  and  F,G  €  M.  The  following  two  state¬ 
ments  are  equivalent: 

1.  dP(F,G)<e 

2.  There  exist  (dependent)  random  variables  U  and  V  with  values  in  fl  such  that  C(U)  =  F 
and  C(V )  =  G,  and  P{d(U ,  V)  >  <}  <  e. 

The  following  results  are  concerned  with  establishing  properties  of  conditional  probabilities 
of  sequences  of  random  variables. 

Lemma  A. 5  Let  Xn  and  Yn  be  sequences  of  random  variables  which,  for  each  n,  are  defined 
on  the  same  probability  space.  Then  Yn  =  op(l)  implies  >  «|Xn}  =  op(l). 

Ptoof  of  Lemma  A. 5:  Y„  =  o,,(l)  implies  that  for  all  «  >  0,  there  exists  N  =  ,N(<)  such  that 
P{||y„||  >  <}  <  fs  whenever  n  >  N.  By  Markov’s  inequality,  for  n  >  N, 

P {p{||y„n  > «|jfB} > e} <  =  flllM > <  f!  =  e. 


Lemma  A. 6  Let  Un,  Vn,  and  Wn  be  sequences  of  random  variables  which,  for  each  n,  are 
defined  on  the  same  probability  space.  If,  for  each  e  >  0,  P{\\Un  —  Vn||  >  e\Wn}  =  or(  1),  then 

dp(£(ffn|vyn),£(yn|vyn))  =  op(i)  • 
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Proof  of  lemma  A.6:  Let  Xn  =  P{\\Un  -  V„j|  >  and  Vn  =  dP(C(UH\Wn),  £(Vn|PPJ). 

Strassen’s  Theorem  implies  that  Yn  <  e  whenever  Xn  <  <•  It  follows  that  P{Xn  >  *}  >  P{Yn  > 
<}  but  P{Xn  >  (}  <  e  for  n  large  enough  and  the  lemma  follows.  ■ 

Lemma  A. 7  If  D^  is  a  positive  random  variable  and  E(Dn\Z„)  =  op(r„),  then  Dn  =  op(r„). 

Proof  of  Lemma  A. 7:  Let  X„  =  E(D„fZ„).  By  definition  ,Y„  =  op(r„)  if,  and  only  if,  for  all 
S,e  >  0  there  exists  N  =  N(S,e)  such  that  for  n  >  N ,  the  set  An  =  {Xn  >  e6rn}  has  probability 
P{An}  <  5.  For  n  >  N ,  we  have  on  A* 

_ _ _  E{Dn\Zn}  ^  Sern  _  , 


P{Dn  >  ern\Zn)  < 


<  -LZZ  =6. 

(  Tn 


Hence, 


P{Dn  >  ern}  =  J  P{Dn>trn\Zn}dP 


<  i  +  8  =  2<5J 


A. 2. 2  Proofs  for  Chapter  3 

Proof  of  Lemma  3.1:  By  Lemma  2.11 


E9%  =  J  9(v)Kn  •  f{u/)dw  =  J  9(u)(f(Jj)  +  Mf  0{\ogn/n))dui 
=  j  6(u>)f(w)dw  +  MfO{\o%n/n)m 


Proof  of  Lemma  3.2:  By  Lemma  2.14 
A„  =  c  ov(9t,9v) 


=  J  J  8(u))d>(u/')cov(I*(u).I%\*'))ditidv' 

=.  ^  J  j9(u)d>(u>')  ^j  F(x)G[T)Dn(x-u>)Dn(x-w')dx  dwdw1 

+— j  J  J  0(w)0(w')  J  J  F(x)G(x)Dn(z  -  v)D„(x  +  w')dz  dudui' 


—  Bn  +  Gi 

Expanding  the  squared  term  in  Bn  above 


Bn  =  J  8{v)  j  J  F(z)G(z)G(y)F(y)Dn(z  -w)Dn(y 

x  f  <t>(<t>')Dn(z  -u')Dn(y-  v')dv' dxdydw . 

J 


By  Lemma  2.10  we  have  that,  uniformly  in  <t>  and  9, 


J  <t>(w')Dn(x  -  u i')Dn(y  -  w')dw'  =  <j>{x)Dn(x  -  y)  +  Af^O^logn), 


and  analogously 


J  9(w)Dn(z  -  w)Dn(y  ~w)cL‘  =  9(z)Dn(z  -  y)  +  MtO( log  n) 
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It  then  follows  that,  uniformly  in  F ,  G,  0,  and  <j>, 

^  f  J  F(z)G(z)F(y)G(y)  (9(z)Dn(z  -  y)  +  MeO(\o g  »))  x 
x  (<p(x)Dn(x  -  y)  +  Af*0(log  n))  dxdy 

-  ~3  f  f  F (x)G(x)F(y)G(y)(<t>(x)e(x)nK*(x  -  y) 

+  G(log  n)(M*0(z)  +  Afg<j>(z))Dn(x  -  y)  +  A/#Af^O(log5  n))  dydx 

=  ~J  Flx)G&jtHzMx) 

+0(logn/n3)  y*  F(*)G(^(Af#fl(*)  +  Af,0(*))  J  J[y)G(y)Dn{x  -  y)dydx 

+0(log3  n/n1)Mg \j  F(x)G(x)dx 

=  1/ /(z)s(x)0(z)^(*)di  +  0(logn/n3)MF^sup  |FG0(/>| 

-t-0(log3  n/n3)  sup  |FG|3  (2  max{Afg  sup  Af*  sup  |6>|}  +  M9M '*) 

Since  #  and  4>  sue  in  Lip(l)  and  are  fixed  functions,  we  shall,  for  notational  convenience,  drop 
uniformity  in  9  and  </>  and  so 

Bn  —  J  f(z)g(x)0(x)4>(z)dx  +  MF^sup  |FGIO(log  n/n3)  4-  sup  |/s|0(log3  n/n3), 

uniformly  in  F  and  G.  The  result  now  follows  immediately  upon  noticing  that  the  bound  on 
Cn  is  the  same.  ■ 

Proof  of  Theorem  3.3:  By  Lemma  3.1  we  have  that 

E\Un~Vn\*  =  nE(9x  -9,  -{8,-8,))2 

=  nVar (9Z  -  9y)  +  n  {e[9x  - 9 ,)  +  E(9y  -  0,)}’ 

=  n  ^Var^,)  +  Var(0y)  —  2cov(0I,  0y)^ 

+(M,  +  Af/)30(log3  n/n) 

so  that  the  theorem  follows  readily  from  Lemma  3.2.  ■ 

Proof  of  Theorem  3.4:  The  regularity  conditions  on  /„  and  Theorem  3.3  imply  E{\Un  ~ 

V„|3|X"  }  =  Op(l)  which,  by  Lemma  A. 7,  in  turn  implies  -  V"n|  =  op(l)  .  ■ 

Proof  of  Corollary  3.5:  (1)  follows  from  the  fact  that  (7„  -  Vn  =  op(l)  implies  (Lemma  A. 5) 

P{\Un-V«\  >  ei-X’o”1}  =  <v0)  which  in  turn  implies  (Lemma  A.6)  dP(C(Un IJfJ-1),  £(V„|A'o  _1)) 
op(l)  but  £(f/n|Xj-1)  =  C(Un).  (2)  follows  from  the  triangle  inequality  and  the  fact  that 
C(Un)  -  C(U).  ■ 

Proof  of  Lemma  3.6:  Hasminskii’s  and  Ibragimov's  (1986)  result  (3.3),  can  be  used  to  estab¬ 
lish 

C{a'y/n(9,  -0/)^  N(0,a'Ea) 

for  any  vector  a  so  that  the  result  follows  by  taking  limits  of  characteristic  functions.  ■ 

Proof  of  Theorem  3.8:  Using  a  Taylor  expansion  we  can  see  that 

;=i  ’ 

=  U"  +  An 


j  F(y)G(y)Kn(z  -  y)dydx 
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where 


An  <  \/n||0*  -  Of  | 


3h(8f) 

3h(8\) 

38 

38 

=  Vn\\Ot  —  Of\\A'n 


and  is  between  0*  and  0/. 
Similarly, 


V-  , 


J=1 


38, 


(8i}  -  Of — )  +  Bn  +  C„ 


=  v;-  +  Bn+cn 


where 


Bn  < 

cn  <  vz\\e9  -  9fn\ | 


II  96 

WfJ 

38 


dh(8D 

38 


9Wl) 

38 


-  \/’»ll Of  -  OfJ\C'n, 


and  0j  is  between  0j  ard  8 j  . 

Now,  Un’  —  V**  =  Op(l)  by  Theorem  3.4  and  by  Theorem  3.7  v/n(S*  —  8f)  has  a  limiting 
distribution  and  hence  ~  Of II  =  Op(l).  Theorem  3.4  implies  that  \/n(8f  -  8 f-  )  has  the 

same  limiting  distribution  as  \/n(S*  —  8f)  and  so  y/n\\8t  —  8 f-  ||  =  Op(l). 

It  remains  to  show  that  A(,,  B'n,  and  C'n  are  each  Op(l).  A'n  =  op(l)  follows  from  ( (0J  ~8f  |[  < 
||0*  -0/H  =  Op(n_1^J)  and  the  fact  that  h  has  continous  derivatives  in  IV#,.  On  the  other  hand, 
if  d(f ,  /„)  =  op(l),  then,  for  n  large,  8 t  £  JV#,  with  high  probability,  and  since  h  has  continous 
derivatives  in  Ngt,  this  implies  that  B'n  =  Op(l)  and  C'n  =  op(l).  ■ 

Proof  of  Corollary  3.9:  Theorem  3.8  establishes  17*  -  Vff  =  Op(  1).  The  rest  of  the  proof  is 
identical  to  that  of  Corollary  3.5.  ■ 


A. 3  Chapter  4 

A. 3.1  Distances  between  spectra 

Definition  A.S  Let  f  and  g  be  two  complex-valued,  measurable  functions  defined  on  [0, 1]. 
Define  the  metric  di(f,g)3  =  /J  |/(u>)  -  g(u/)\3du>. 

Lemma  A.S  If  f  and  g  are  complex-valued,  measurable  functions  on  [0, 1],  then 
dj(f,g)  <  sup  max(  (/(<*>  )|,  |j(in)|)dj(log  /,  log  y). 

(a/ 

The  proof  is  a  consequence  of  the  following  lemma: 

Lemma  A. 9  |1  —  e*|  <  |z|max(|e*|,  1)  for  complex  z. 

Proof:  Let  z  =  x  4-  iy.  Then 

\e‘  -  1|J  =  eu  —  2e*  cos(y)  +  1 

=  (e*  -  1)J  +  2e*(l  -  cos(y)) 

<  (e1-!)3 +eV 


(A.6) 

(A.7) 

(A.S) 
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since  1  -  cos(y)  <  y3/2.  By  Taylor’s  theorem,  there  exists  0  <  9  <  1  such  that  (e*  -  l)3  = 
z3e#3*  <  x3max(e3*,  1).  Hence  |e‘  —  lj3  <  (x3  +  y3)  max(e3*,  1).  The  result  follows  upon  taking 
square  root  and  using  the  fact  that  ]e*|  =  ez .  ■ 

We  can  now  prove  lemma  A. 8:  Apply  Lemma  A. 9  with  z  =  log  g  -  log  f: 

\f  -  9l  =  1/IU  -  ?// 1  <  l  log  /  —  log  s|  max(|/|,  |s|).  (A. 9) 

The  result  now  follows  upon  integrating  A. 9.  ■ 


A. 3. 2  Proofs  for  Chapter  4 

Proof  of  lemma  4.1:  Lemmas  2.11  and  2.13  show  that 


Vkn=(l+Pk0k)(U?k+Utk)/2 

where  Ulk  are  iid  N(0, 1)  variates,  6k  is  a  random  variable  such  that  |0*|  <  1,  and  pk  is  the 
correlation  between  the  real  and  imaginary  parts  of  the  DFT.  Thus  we  have  that 

(kn  —  log(t/3*  +  U\k)/2  4-  log(  1  +  Pk&k)  +  7- 

Lemma  2.12  establishes  pk  =  M/0(log  n/n).  Since  (U*k  +  U%k)/2  is  an  exponential  random 
variable,  the  expected  value  of  its  log  is  just  — y  so  that 

Etkn  =  pkE0k  -  p\E6\l2  +  p\E0\! 3  +  0(p4k). 

The  first  assertion  then  follows  from  the  fac»  that  E(&™)  —  0  for  m  odd  and  the  second  follows 
trivially. 

Now,  the  variance  of 


=  l°g(l  +  Pk#k)  +  l<>g(iu  +  U\k)/2  +  1 


is  equal  to  the  sum  of  the  variances  of  the  first  two  terms  since  lemma  2.13  shows  they  are 
independent.  The  second  term  is  just  the  log  of  an  exponential  so  its  variance  is  known  to  be 
t3/6.  Moreover,  the  variance  of  the  first  term  is  bounded  by  (log(  1  -  pk))1  =  0(p\)  so  that  the 
third  assertion  follows  trivially. 

The  fourth  assertion  is  merely  a  restatement  of  lemma  2.17.  ■ 

Proof  of  lemma  4.2:  Notice  that 


SUP  $ 


2 -6k  _  (2ra)»  2 

»  »>o  1  +  A(2 ri)3™  +  n  J0 


(2  ■Kx)kdx 
1  +  A(2rr*)3m 


The  result  now  follows  upon  changing  variable  of  integration  to  y  =  Al/3m2i rx  in  the  integral 
above  and  verifying  that 


®>S  1  +  A(2jtz)3"*  \A(2m -Jfc), 

Proof  of  Lemma  4.3:  By  Lemma  4.2  and  since  /  €  Lip(l)  we  have  that,  for  n  large, 


(2**)k 


(k  \  kt,m 
M2^T))  =0(a 


-fc/3m 


n 

A»M  =  ^2  (log/(*/2n)  +  O(logn/n)}  W(u)  -  i/2n) 

*=-n+l 

n 

=  J2  log/(ife/2n)W(^  -  fc/2n)  +  (l +A-l/3m)0(logn/n) 

*=-n+l 

=  £«(«)  +  (!  +  A_l/3m)0(log  n/n), 
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so  that  it  suffices  to  show  supw  |Bn(<»))|  <  |fc„|.  Indeed, 


W  =  2^  £  ^  exP(»2irt(i/  —  s)/2„) 

i  =  -n  '  '  fc  =  -n-f  1 


En  6„exp(«2Ti/w) 
1  +  A(2*i/)lm 

—  n  v  ' 

oo 

<  £  IM-« 


(A. 10) 


Proof  of  Lemma  4.4:  A  similar  derivation  to  that  used  in  deriving  (A. 10)  in  the  proof  of 
Lemma  4.3,  shows  that 

jL(w)=  V  ^n«p(,_2j^) 

1  1  A  1  +  A(2iri/)Jm 

where  {6„n}  the  Fourier  coefficients  of  log(JfJn  *  /).  This  implies  that 

KM  =  2nW\,m,n  .  log (Kin  .  /)(<*) 
and  hence  by  Lemmas  2.5,  2.6,  and  4.2,  we  have  that 

sup  KM  I  <  2nsup  KJn*  ||V^X.m.nlll 

<  sup|/'(u>)|sup  (1  +  A*1/3m)0(l). 

w  w  /M 

uniformly  in  A.  But  /'  is  bounded  since  /  €  Lip(l)  is  differentiable.  Furthermore,  since  /  >  0, 
and  is  continuous,  1//  is  bounded.  ■ 

Proof  of  Lemma  4.5:  The  r.h.s.  in  4.6  follows  from  the  fact  that  (cf.,  e.g.,  the  proof  of  Lemma 

4.3) 

*•(«)  =  L  +  t1  +  A ~ 1/3m )0(log  n/n), 

Parseval’s  theorem,  and  the  fact  that  (cf.,  e.g.,  Butzer  and  Nessel,  1971) 


£  =  llfljn+i  •  A  -  h||j  =  0(n-s) 


whenever  h  €  Lip(l).  Finally,  continuity  of  A(A)  follows  from  a  discrete  version  of  Lebesgue’s 
dominated  convergence  theorem.  ■ 

Proof  of  Lemma  4.6:  From  (4.4)  and  (4.5)  we  have  that 

1  V'  V"  exP(i2irt/(<x  ~  */2n)) 

kn^’  ~  2n  ^  *n  2-  1  +  A(2ri/)3m 

h= -n-f  1  v=-n 

_  A  exp(»2TiM  /  1  A  _ , 


=  V  j  _L  V  fkn  exp(— »2»i/fc/2r 

^  1  +  A (2xu  ,m  \  2n  t  ^  FV  ' 

t/=  ~n  '  '  \  fc  =  -«  + 1 

E"  exp(»2irvw)B  lstl 
1  +  A(2M,m 
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Now,  by  Lemma  4.1 


£!B„n|J  =  E  E(fknejn)exp(-i2iru{k  -  j)/2n) 

],k=-n+l 

<  0(n~l)  4- 0(log3  n/«3)  =  0(n-1) 


uniformly  in  u. 

On  the  other  hand,  by  Parseval’s  Theorem 


/  E{h-n(«?}cL,  =  £  ~ 

*  !/  =  —  »!  ' 


=  E  7T 


=  0(n'l)(l  +  A-1/3m) 

by  bounding  the  sum  by  an  integral  and  changing  the  variable  of  integration  as  in  the  proof  of 
Lemma  4.2.  ft 

Proof  of  Theorem  4.7:  Take  0  <  »  <  1  and  divide  up  the  2n  integers  from  — n  +  1  to  n  into 
J  ~  n1"  consecutive  disjoint  sets  E ,  from  ai  +  1  =  -n  4-  1  to  aj,  aj  4-  1  to  aj,  etc.,  each  set  of 
length  n‘  or  less.  Let 

=  ^  1  etn, 


1  <  j  <  J ■  Then, 


.a* -o 


Var(i4,„)  <  n'O(l)  4-  n3'0(log3  n/n3)  =  O(n'), 
vhere  0  is  uniform  in  j  for  1  <  j  <  J.  Applying  Markov’s  inequality  to  Aj3n  we  have  that 


e{u,„|  >  .■/■**}  <  =  o, 

and  if  An  =  maxi <;<r  |A>n|  then 

F{An  >n1/3+1'}  <  £>{[Ajn|  >  n1/3+<} 


_  W)  QK-v)  u 

~  ni+j«  nl+3<  ~ 


Thus  An  =  of,(n1'/3+<),  for  each  S  >  0.  Now, 

j 

KW  =  E  E  W*(U  -  */2n)‘*" 

=  £  J  W»(w  ~  °;/2nM;»  +  E  “  */2n)  -  ~  «j/2 »))e*n 

r=t  I  MS, 


appendix  a.  auxiliary  results  and  proofs 


The  first  sum  above,  Si„,  satisfies 

J 

Sin  <  A„sup^  |WB(w  -a,/2n)\ 

"  3=1 

<  A„  sup  |Wrn(<ii)|  J 

W 

=  op(nl/i+‘)Olnr'-l)0(n1-) 

=  Op(n1,1+‘+r°-‘) 

and  the  second  sum,  S2n,  satisfies 

j 

Sin  <  ~  fc/2n)  -  W’nO*'  _  a;/2n)Hfknl 

j=l  k €E, 

<  ^2  — -0(n'‘~l)  5^  l«*»l 

J  =  l  ”  k€Bj 

<  -o(«'‘-1)E  £  i“»i 

n  3  tee, 

=  nOp(n‘+,‘ _J) 

=  o,(«*+'‘-1) 

since  Y.k  I4*"!  ~  °i>(n)  by  Markov’s  inequality.  This  proves  the  first  assertion  of  the  lemma. 
An  analogous  argument  proves  the  second.  M 

Proof  of  Lemma  4.8:  By  the  triangle  inequality 

di(h,hn)  =  di(hX+K) 

<  di(h,hn)Adi(0,h-J 

and  the  result  follows  immediately  from  Lemmas  4.5  and  4.6  since  Edi(h’n,  0)  < 

-/EdffaTo).  .  .  ■ 

Proof  of  Lemma  4.9:  Apply  Lemma  A. 8  to  /  and  /„.  The  Or  term  follows  since  /«  is 
bounded  in  probability  (Lemma  4.7).  The  expression  in  parentheses  follows  from  Lemma  4.8.  ■ 

A. 4  Chapter  7 

Proof  of  Lemma  7.1:  Since  V  =  mod(f7  +  a,  1)  ~  U(0,  1)  for  all  a  whenever  U  ~  17(0, 1),  we 

have  that  _ 

y"-1  =  Fnti'^V 

haa  the  same  distribution  as 

FX/JDBr;-1 

where  Z)n  is  a  diagonal  matrix  with  entries  exp(i2wt7o),  •  •  • ,  exp(»2x(7n_,),  U'i  iid  17(0, 1)  so  let 

Ko"-1  =  KAn^D*^-1 

=  FnDnA'J’Fo-1 

=  FnDnFnX^-1 


so  that  the  covariance  matrix  of  '0"_1  given  Dn  is  given  by 
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On  the  otner  hand,  for  any  matrix  A«  =  (a,*),  DnAnDn  =  ( a,k  exp(»2x(t/j  -  17*)))  .  and  s‘nce 
E txp(i2x(U]  -  Ui,))  =  1  only  if  j  =  k  and  0  otherwise,  EDnAn Dn  —  diag{An),  i.e.,  only  the 
diagonal  survives  premultiplying  by  Dn,  postmultiplying  by  Dn,  and  taking  expectation.  Thus, 
taking  expectation  again  in  (A. 11)  yields 

£{yo"-1(yon-,n  =  T;diag(FnE„T;)F„ 

=  e; 


where  is  as  in  Example  1  ofcirculants  in  Chapter  2. 
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estimate  that  under  /.  We  establish  results  for  estimators  of  well-behaved  functions  of  linear 
functionals  of  the  spectral  density.  In  the  second  problem  we  pretend  the  data  are  circular  so 
that  the  periodogram  of  the  data  is  a  sufficient  statistic  and  thus  the  conditional  expectation  of 
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